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Numerieke approximatie. 


Bij tal van mathematische problemen is de exacte oplossing niet 
in handzame vorm te vinden. Zo is fEGOdx maar voor relatief 
weinig functies f als getalwaarde be bepalen, zijn differentiaal- 
vergelijkingen nog zeldzamer in eindige termen van elementaire 
functies oplosbaar en zijn wortels van vergelijkingen f(x) = 0 
slechts zelden te bepalen. En de oplossingen van deze problemen 
worden zelfs ongedefiniëerd als men van de optredende functies 
slechts in een aantal punten de waarden kent, terwijl dit bijv. 


bij de experimentele wetenschappen toch een veel voorkomende 


situatie is. 


Asymptotische ontwikkelingen en de Ritz-Galerkin methode zijn 
voorbeelden van methoden om benaderde uitkomsten voor de twee 
eerstgenoemde problemen te vinden. In het huidige college zullen 
wij laten zien hoe men de genoemde problemen benaderd kan oplossen 
door de optredende functies te benaderen door functies uit 

zinnige functieklassen waarvoor de gestelde problemen wèl op- 


losbaar zijn. 


Wanneer men een functie f zo goed mogelijk wil benaderen door een 
functie uit een zekere functieklasse dan moet men op een of andere 
wijze een afstandsbegrip in de functieruimte definiëren (bijv. 

als norm van de verschilfunctie) en dan trachtendie g te vinden 
waarvoor afst(f;g) zo klein mogelijk is. Dit pleegt nogal een 


werk te zijn. 


Wanneer men evenwel een functie slechts op een klein interval wil 


approximeren doet men vaak goede zaken met polynoom interpolatie, 


d.w.z dat men een polynoom zoekt dat in een aantal gegeven punten 


met de te approximeren functie overeenstemt (wat van het zojuist 
genoemde approximatie standpunt uit bezien eigenlijk een absurde 
eis is), ook al vanwege de gemakkelijke construeerbaarheid van 
deze approximatie. In het pre-computer tijdperk bestonden die 
kleine intervallen vaak uit stukjes tabel, en benaderde men m.b.v 
interpolatie tussengelegen waarden. Dit is nu wel enigermate uit 
de tijd, maar benadering op kleine intervallen is nog steeds 

van belang om de in de aanhef genoemde problemen benaderd op te 
lossen. Vandaar nu eerst een hoofstuk waarin interpolatie nader 


bestudeerd wordt. 
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A. 


Hoofdstuk 1. Interpolatie. 


Inleiding 

Het berekenen, of zelfs het enigszins redelijk benaderen, van 
de waarde van een funktie in een punt is vaak erg moeilijk of 
tijdrovend; zeker als er veel funktiewaarden worden gevraagd. 
Het ligt dan ook min of meer voor de hand dat men voor slechts 
enkele punten de funktiewaarde bepaalt en daaruit op een goed- 
kope manier een benadering van de funktiewaarden in de andere 
punten probeertaf te leiden. Bij interpolatie komt deze goed- 
kope manier erop neer dat men de funktie benadert door een 
polynoom (een zgn. interpolatiepolynoom) dat in een aantal 
punten (de steunpunten) dezelfde waarde aanneemt als de ge= 
geven funktie. We illustreren dit idee aan de hand van lineaïre 


interpolatie. 


Lineaire interpolatie 


Laat a en b twee punten zijn in IR, zeg a <b, en f een funktie 
die gedefinieerd is op [a,b]. Het lineaire interpolatiepolynoom 


van f met steunpunten a en b wordt dan gegeven door 
ks) px) = fla) + FED) - f(a). 


Grafisch komt dit erop neer dat men de grafiek van de funktie 


benadert door de rechte lijn door de punten (a,f(a)) en (b,f(b)). 
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Lineaire interpolatie is een veel nauwkeuriger proces dan de 
meeste mensen denken. Beschouw bijvoorbeeld nevenstaand tabelletje 


van de sinus. Stel eens dat we sin(37°) 


niet gekend hadden en deze door lineaire X sin(x) 

interpolatie uit sin(36° ) en sin(38°) 35° 0573576 
hadden willen bepalen. Dan was er uit- 36° 0.587785 
gekomen 0.601723, hetgeen slechts 9 een- gj 0.601815 
heden van de vijfde decimaal verschilt van 38° 0.615661 
het ware antwoord. 39° 0.629320 


De interpolatiefout loopt snel op bij het verlengen van het 
interpolatieinterval (en daalt dus snel bij het verkorten!). 

Zo heeft sin 37° , berekend m.b.v. lineaire interpolatie op sin 35° 
en sin 39° al een fout van 37.107". Dat de fout bij verdubbeling 
van het interpolatieinterval ongeveer vier keer zo groot is ge- 
worden, is geen toeval. Om dit in te zien leiden we nu een uit- 


drukking af voor de fout. 


Eerst schrijven we (1.1) om tot 
_ b=-Xx x-a 
(1.2) plx) = 5ozfla) + Sgf(b) 
Als f een continue tweede afgeleide heeft op [a,bl, en x tussen 
a en b ligt, volgt uit de formule van Taylor 
_v2 
Fla) = EG) + (amd!) + BEL ere), a SE, Cx 


en 
(bx)? 


F(x) + (b-x)f!'(x) + 5 


f(b) 


FRLE De XS Er SB 


Invullen in (1.2) geeft 


5 _ Lbexilarxt (Ra oa DX en 
px) = f(x) 5 [5-5 FCE, + pit (Bid 


Nu staat er tussen [...l een gewogen gemiddelde van £"(E,) 


en f"(E,) = dat is dus een getal tussen f"(E,) en f"(E,) -— 


1. 


en wegens de doorlopendheid van f" bestaat er dus een £ tussen 


a en b zodat 


u _ XTA en b-Xx “ 
LUIE) = Bed f (E‚) + DE 5 CE) 


We vinden zodoende 


Stelling (1.3). Als f op [a,bl een continue tweede afgeleide 
heeft en p(x) het lineaire interpolatiepolynoom met steunpunten 


a en b is, dan bestaat er bij iedere x Ela,bl een £ Ela,b) zodat 


EGO =p) = Graden) LE 


In het voorbeeld, waar we sin 37° benaderden m.b.v. lineaire 
interpolatie op 36° en 38° is f"(E) = -sin E en E ligt tussen 
36° en 38°. Dus is in (1.3) f"(E) = -0.6. 

Verder is x-a = b-x = 1° = El radiaal. De formule in (1.3) 

geeft nu f(x) -— p(x) # 0.000092, hetgeen precies klopt met de 
gevonden fout van 9 eenheden in de vijfde decimaal. 

Merk op dat het nu ook duidelijk is, waarom de fout bij verdubbe- 
ling van het interval vier keer zo groot werd. 

De absolute waarde van de uitdrukking (x-a)(x-b) is, wanneer we 
alleen x-waarden tussen a en b bekijken, hoogstens gelijk aan 
A(b-a)® 5 dit maximum wordt aangenomen voor Xx = (atb). 

Derhalve 

Stelling (1.4). Bij lineaire interpolatie is de fout in abso- 
lute waarde hooguit 


5b-a)? max TELE) 
EE[ a ‚bl 


Interpolatie in het algemeen. 


Algemenerstelt men zich bij interpolatie de opgave bij een 





gegeven funktie f en gegeven (steun)punten XgsrrrsX, een 

hoogstens n-de graads veelterm p te bepalen zodat f en p in de 
gegeven punten overeenstemmen. 

Alhoewel we er in het volgende slechts incidenteel gebruik van 
zullen maken bekijken we nu een nog wat algemenere situatie nl. dat 
ook nog afgeleiden van f en p in de gegeven punten overeenstemmen. 
Laat Xpo Xn verschillende punten zijn en Pose: or, een stel 


niet-negatieve gehele getallen. Zij f een funktie die in X; minstens 


r;‚-maal differentieerbaar is (Ì = 0Os...‚n). Een polynoom p met de 
eigenschap p (xj) = Eed voor k = 0,...r; en i == Os...‚n heet 
dan een interpolatiepolynoom vän f met steunpunten Xs... ;X« 


n 
Het is de bedoeling de graad van p zo laag mogelijk te houden. 
Welnu, de voorwaarden dn, = ges leveren N + 1 lineaire 
vergelijkingen voor de coëfficiënten van p (vgl.(2.2)), waarbij 
N= nt rg + eee + Et We proberen daarom eens een polynoom 


met N + 1 coëfficiënten, d.w.z. een hoogstens N-de graads poly- 


noom, te vinden. De volgende stelling laat zien met welk resultaat. 


Stelling (2.1). In de hierboven beschreven situatie bestaat er 


precies één interpolatiepolynoom p waarvan de graad < N is. 


Bewijs. Zij plx) = ae & li +... t a. De eisen 
pd = fra) voor k = Oss sr; en Ì = 0,...‚n geven 


aanleiding tot N+1 lineaire vergelijkingen met de coëfficiënten 


Ag». >äN als onbekenden bijv. 
Ned n 
An*i t Ay-4 Xi Heee + a, = f(x) 
€22) 
N= 1 ie N-2 En r > 
Nayx; + (N-Day4*; tee ta, = f'(x) als r; > 1 


De bewering van de stelling is, dat dit stelsel precies één op- 


lossing heeft. De matrix van dit stelsel hangt niet af van de 


waarden van f en eventuele afgeleiden van f. We kunnen daarom 
volstaan met aan te tonen dat het homogene probleem, d.w.z. 


f == O0, slechts de nuloplossing heeft. 
(k) 


De eis p (xj) = 0 voor k = 0,...sr; impliceert dat p deelbaar 
„+1 
is door (x-x,) Ì . Dit geldt voor i = 0,...‚n en dus is p deelbaar 
ry +1 Pati 
door Kek)  auvaas rx) . Deze faktor heeft evenwel graad 


N+1, terwijl de graad van p hoogstens N is. Derhalve kan p hoog- 
uit het nulpolynoom zijn. Het is onmiddelijk duidelijk dat het 


nuipolynoom inderdaad ook voldoet. 


Opmerkingen. Als we in het volgende spreken over het interpolatie- 
polynoom, dan bedoelen we dat met graad SN. 

Merk op dat het aantal coëfficiënten van p hoogstens gelijk is 

aan het aantal eisen waaraan moet worden voldaan. Met dit tel- 
argument kan men de stelling gemakkelijk onthouden. Het inter- 
polatiepolynoom kan best graad < N hebben. Dat is bijv. het geval 
bij het homogene probleem hierboven. 


We geven nu een uitdrukking voor de fout bij deze algemene inter- 


polatie. 
Stelling (2.3) Zij la,bl een segment dat X,s:.:»X, bevat en laat 


p het interpolatiepolynoom zijn. 


(N+1) 


Veronderstel dat f bestaat op (a,b). Dan is er voor elke 


x E{a,bl een E E(a,;b) zodat 


Pr td Pyt roti N41, ) 
Fl) pix = Gad . (ex) eene rx) EEE 
Bewijs. Als x samenvalt met een van de punten Xpo esXn kan 


men & willekeurig kiezen. 





Laat nu x # X; zijn voor alle i. In de rest van het bewijs is x 
vast ; t is de variabele. 
Bekijk de funktie & met 


r.+1 


El ‚+ Ì 


n EL Ë 
BE) = CFE) — plt)) MH (x-x‚) > Gfla) > plaid TL (Em) E 
i=0 i=0 


Het is duidelijk dat KX ee Xr nulpunten zijn van &; dus 

heeft & minstens n+2 nulpunten. Met Rolle volgt hieruit dat $'!' 
minstens n+1 van Xx en de Xj verschillende nulpunten heeft. Boven- 
dien heeft 6 nog een hulu in x als r;>1 is; zeg dat r, > 1 
optreedt voor k waarden van i, dan geeft dit nog eens k nulpunten. 
In het totaal heeft &' dus minstens nt+k+1 nulpunten. Er volgt, 
wederom met Rolle, dat &' minstens n+k nulpunten heeft verschillend 
van Xgsee-»Xjr Daarbij komt nog een nulpunt in iedere x; waarvoor 


(N+1) 


p, > 2 is. Zo doorgaand kan men inzien dat ® minstens één 


nulpunt & heeft op (a,b). Deze E is de gezochte ; immers, omdat 


n Pet 
de (N+1)ste afgeleide van p nul en die van Il (t-x;) ú gelijk 
i=0 
aan (N+1)! is hebben we 
n r.+1 
vat gps UE T (ex) È (EGO =- PGO NH)! 


1=0 


Opmerkingen. a. Door in het voorgaande n=0 te nemen, vindt men 
als interpolatiepolynoom het begin van de Taylorreeks van f (tot 
en met de r‚-de graadstermj 5; de formule voor de fout, die in 
(2.3) is afgeleid, geeft dan de bekende restterm van Taylor. 


b. Wanneer men in het voorgaande Lg °° By Faes Ei B 0 neemt, krijgt 


men zgn. Lagrangeinterpolatie ; men eist dus slechts dat het 


interpolatiepolynoom in de punten x > Xn dezelfde funktie- 


02 


waarden heeft als f. Lineaire interpolatie is hiervan een 


KE Oelimenie 
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bijzonder geval. De formule voor de fout zoals die in (1.3) 
wordt gegeven,is hetzelfde als die in (2.3) ; merk echter op 

dat in (2.3) x niet persé in het door de steunpunten opgespannen 
segment hoeft te liggen. 


ce, Wanneer men Pr, = Pj Zes. = P_= 1 neemt, krijgt men zgn. 
0 1 n 


Hermiteinterpolatie 


Lagrangeinterpolatie. 


We adviseren niet om het interpolaätiepolynoom ook inderdaad te 
bepalen door oplossing van het lineaire stelsel dat in het be- 

wijs van stelling (2.1) optrad. Dat is namelijk tijdrovend, terwijl 
oplossing van het betreffende stelsel erg gevoelig voor afrond- 
fouten pleegt te zijn. 

Er bestaan allerlei andere manieren om het interpolatiepolynoom 

te bepalen (zie bijv. 84). Een manier waarvan men nog wel eens 
plezier kan hebben in het belangrijke geval van Lagrangeinter= 


polatie is een generalisatie van (1.2) 


& (n) 
erk plx) = Ge Le Cad EC) 


k 


waarbij 


(BaZ) 





(de (n) duidt hier geen afgeleide aan!) 
De juistheid van (3.1) volgt onmiddelijk uit de opmerking dat 


ed = 0 als Ìì £ k en tia = 1. Men noemt deze schrijf- 


wijze de Lagrangerepresentatie van het interpolatiepolynoom; de polynomen 


Aad heten de Lagrangecoêfficienten. 


Een belangrijk geval krijgt men als de afstand tussen twee op- 
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eenvolgende steunpunten steeds gelijk is, zeg gelijk aan nh. Dan 


nummert men de steunpunten zo dat Xi 5 Xo + kh voor k = 0s...jn. 
Men noemt dit equidistante Lagrangeinterpolatie. Schrijven we 
nu = X4 + ph, waarbij i een van de indices 0,...‚n is, dan 


krijgt men 


(3.3) AND. 2 pri) prionen (ptiekri) prik) BEN (peeler) 
Kel) ki (n=k)S 
voor k = 0O,...‚n. Dit is blijkbaar onafhankelijk van h. 
Voor dit geval zijn tabellen van de Lagrangecoëfficienten voor- 
handen, hetgeen plezierig is voor het rekenen met de hand of 
met een tafelrekenmachine : zie bijv. voor n=3 : Abramowitz-Stegun, 
Handbook of Mathematical Functions, p901. Deze tabel hebben we 
op de volgende bladzijde overgenomen. Hiermee krijgt men bijv. 
sin(37.01° ) = -0.0832835 sin 36° + 0.9949005 sin 37° + 0.0100495 
sin 38° — 0.0016665 sin 39° 
(neem i = 1 en p= 0.001) 
en. sin(37.99°) = -0.0016665 sin 36 + 0.0100495 sin 37° + 0.9949005 
sin 38° — 0.0032835 sin 39° 
(neem i = 1 en p= 0.99). 
Volgens (2.3) bedraagt de fout bij (n+1) -punts (oftewel n-de orde) 


Lagrangeinterpolatie 


(nt op) 
(3.4) CRH NRE Ae ane et (xr x) CRETE 


Bij funkties die men in de praktijk tegenkomt zijn de hogere af- 
geleiden vaak bijzonder ingewikkelde uitdrukkingen, waarbij het 
vrijwel ondoenlijk is te schatten hoe groot de waarden wel kunnen 
worden. We zullen overigens aan het eind van deze paragraaf en ook 
in de volgende zien hoe men zieh toch nog een indruk kan verschaffen 


van de waarde van et), 


he 


NUMERICAL ANALYSIS 


FOUR-POINT LAGRANGIAN INFERPOLATION COEFFICIËNTS 


A 


0. 00000 
-0, 00328 
-0, 00646 
0. 00955 
-0. 01254 


-0. 01543 
-0, 01823 
-0, 02094 
-9. 02355 
-0, 02607 


-0, 02850 
-0, 03083 
-0, 03308 
-0, 03524 
-0, 03732 


0, 03931 
-9, 04121 
-0, 04303 
0, 04477 
-0, 04642 


-0, 04800 
-0, 04949 
-0, 05090 
-0, 05224 
„0, 05350 


-0. 05468 
=0, 05579 
-0, 05683 
=0, 05779 
„0, 05868 


0, 05950 
-0, 06024 
-0, 06092 
-0, 06153 
-0. 06208 


0, 06256 
-0, 06297 
-0, 06332 
-0, 06361 
-0, 06383 


-0, 06400 
-0, 06410 
-0, 06414 
„0, 06413 
-0. 06406 


-0, 06393 
-0. 06375 


-0, 06352. 


-0, 06323 
-0. 06289 


-0, 06250 
Âz 


Áo 
1. 00000 
0, 99490 
0. 98960 
0. 98411 
0. 97843 


0. 97256 
0. 96650 
0. 96027 
0. 95385 
0. 94726 


0. 94050 
0. 93356 
0. 92646 
0, 91919 
0. 91177 


0. 90418 
0. 89644 
0, 88855 
0. 88051 
0. 87232 


0. 86400 
0, 85553 
0. 84692 
0, 83818 
0, 82931 


0, 82031 
0. 81118 
0, 80194 
0, 79257 
0, 78309 


0. 77350 
0. 76379 
0. 75398 
0. 74406 
0. 73405 


0. 72393 
0. 71372 
0. 70342 
0. 69303 
0, 68255 


0, 67200 
0. 66136 
0, 65064 


0. 63985 


0, 62899 


0, 61806 
0. 60706 
0, 59601 
0. 58489 
0. 57372 


0. 56250 
A; 


Afps=t- 12 


00 
05 
40 
35 
20 


25 
80 
15 
60 
45 


A1 
0, 00000 
0, 01004 
0. 02019 
0. 03043 
0. 04076 


0. 05118 
0. 06169 
0. 07227 
0. 08294 
0. 09368 


0. 10450 
011538 
0. 12633 
0. 13735 
0. 14842 


0. 15956 
0. 17075 
0. 18199 
0. 19328 
0. 20462 


0, 21600 
0. 22741 
0. 23887 
0. 25036 
0. 26188 


0, 27343 
0. 28501 
0. 29660 
0. 30822 
0. 31985 


0. 33150 
0. 34315 
0. 35481 
0. 36648 
0. 37814 


0. 38981 
0, 40147 
0, 41312 
0. 42476 
0. 43639 


0. 44800 
0. 45958 
0. 47115 
0. 48269 
0, 49420 


0. 50568 
0. 51713 
0. 52853 
0. 53990 
0. 55122 


0. 56250 


Ao 


pip: (Pp 2) 
CATZ — be) Np —k) 


00 
95 
60 
65 
80 


A2 

0. 00000 
-0. 00166 
-0, 00333 
—-0, 00499 
-0. 00665 


-0, 00831 
-0, 00996 
-0. 01160 
-0. 01324 
-0. 01487 


-0. 01650 
—0, 01811 
-0. 01971 
-0, 02130 
-0, 02287 


-0. 02443 
-0. 02598 
-0, 02751 
-0, 02902 
„0. 03052 


-0, 03200 
-0, 03345 
„0. 03489 
-0, 03630 
„0, 03769 


„0, 03906 
-0, 04040 
0. 04171 
-0. 04300 
-0. 04426 


—0. 04550 
-0. 94670 
0, 04787 
-0. 04901 
-0. 05011 


-0, 05118 
-0, 05222 
„0, 05322 
-0. 05418 
-0. 05511 


-0. 05600 
-0, 05684 
„0, 05765 
-0. 05841 
-0. 05913 


-0. 05981 
-0. 06044 
-0. 06102 
„0. 06156 
„0, 06205 


-0, 06250 
û} 
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Anderzijds; wanneer gint) 


slechts weinig varieert, hetgeen geen 
onredelijks aanname is aangezien men meestal interpoleert op kleine 


intervallen, wordt het gedrag van de fout op het interval bepaald 


door het polynoom rx dese err). We bekijken dit gedrag eens 
in het equidistante geval. Substitutie van x = Xx, t+ ph geeft 
Gez) rx de erx) = nP*ip(p-1).... (pen). 

Voor n = 1,2,3,4,5 zijn de grafieken van p(p-1)....(p-n) hier- 


onder geschetst. 








di 
E 
Or 5 € 
O2 
k 5 
ad 
Ale & 
_ÎÌ 


\N 


AN DV AC NO 








De relatieve extrema worden telkens aangenomen in de buurt van 


de punten 55 etc. De waarden van de extrema zijn er in de 


3.4 5,6 


tekening bijgeschreven. Het blijkt dat de extrema, en dus ook 

de interpolatiefouten, het kleinst zijn midden op het inter- 

polatieinterval; iets waarmee men bij de keuze van de steunpunten 

rekening moet proberen te houden. 

We duidden voorheen al aan dat men zieh toch een indruk kan ver- 
(n+1) 


schaffen over de waarden van f . Een van de methoden daarvoor 


werkt als volgt 


kn har edn a ae 


Stel : We kennen de funktiewaarden in een (groot) aantal punten 


Bn neske reen ‚ (die we eenvoudigheidshalve in volgorde 


oer Knee? 
van grootte genummerd denken : Xx, Hi Swane de Laat Pp, het n-de 
orde interpolatiepolynoom zijn op XporerXjD dan kunnen we de 


waarde van het polynoom en die van de funktie vergelijken in het 


punt Vd” We vinden zodoende de waarde van g{n+i) in ‘rn 
punt Es op (xps, 41 
(n+1) Crrije Nj 
GEEN f CE) = (EME LEM EETL 0) [ECx 44) Po naa ’] 


n+1 


Evenzo vindt men uitgaande van het interpolatiepolynoom p‚ op 


XjorresXn4g CN de funktiewaarde in X42 de waarde van Mt) in 
1 

n &, op Xx, oXn42] . etc. 

Stel nu eens dat deze aldus verkregen waarden de een 


gladjes verlopende rij vormen en laten we deze, althans in gedachten, 
eens grafisch uitzetten d.m.v. horizontale strepen ter lengte 


n+1 op de hoogte de ) boven het interval [ x, ieXzar) Zie 
fig.£.6 voor n=2) . Dan snijdt de grafiek van Inti) an ze 







ECE) 


EEECE 


ÊiE Aas 








-_ 14 =— 


horizontale lijnstukjes en men verwacht dat hij binnen de aan- 


gegeven strook loopt. Op die manier krijgt men dan een indruk 


+ 
van de waarden van ‚mn DD, 


Wie dit een bedenkelijke gang van zaken vindt heeft gelijk, 
maar hij bedenke tevens dat men hier in feite hetzelfde doet 
als wanneer men, om enig zicht op het verloop van een lastige 


functie te krijgen, een aantal funktiewaarden bepaalt (zie fig. 


87) 





„Fina Âa] 


en dan aanneemt dat het verloop door de getrokken kromme wordt 


weergegeven en niet door de gestippelde. 


Gedeelde differenties 


De methode met de Lagrangecoëfficienten is, althans in het 
equidistante geval, handig als men een tabel met Lagrange- 
coëfficienten en een tafelrekenmachine ter beschikking heeft. 

De! methode is ook geschikt voor gebruik op de computer; aangezien 
de Lagrangecoëfficiënten nogal efficiënt te berekenen zijn. 


(Het opnemen van een tabel in het computergeheugen zou te veel 


‚ plaats innemen). Een bezwaar is echter dat men van tevoren de 


orde moet vinden die een voldoend kleine fout geeft. Daarbij 
moeten hogere afgeleiden worden geschat, hetgeen meestal on- 
doenlijk is. Wellicht denkt men daar $3C voor te kunnen ge- 
bruiken, maar dat is veel te omslachtig. Voor ieder lijnstuk 
in fig.3.6 moet men immers een interpolatiepolynoom berekenen. 


Het zou ideaal zijn als men een procédé had waarbij men een- 


voudig de orde zou kunnen verhogen en daarbij gaandeweg een in- 
druk zou kunnen krijgen van de fout; men kan dan stoppen als de 
geschatte foutklein genoeg is. Er zijn inderdaad een aantal van 
dergelijke methodes bedacht; we behandelen hier de methode van 
de gedeelde differenties. 

Zij XgsXjs--… een rij punten die we als steunpunten kunnen ge- 
bruiken; het is niet nodig te eisen dat dit een equidistante rij 
is, of dat de x'en zijn genummerd in volgorde van grootte. Zij 


(ij het interpolatiepolynoom op Xoor ersXj-g EN Ë dat op Xx 


n=1 en 


02e 
Dan is Ea Ere een hoogstens n-de graads polynoom met Xp orerrkned 


als nulpunten. Dus is Ea ed een constante maal (Xx) xp es. rx) « 


Deze constante noteren we als fXgse rex) ; het is blijkbaar 
de coëfficiënt van Xx in fx). 


Het is duidelijk dat 
(4.1) ErOodsf dt erg df Og ox tang dax FCX ox oa dte eet 
trg dees erg Eg sere sj) 


We bekijken nu hoe we hier in de praktijk gebruik van kunnen 
maken. 


Schrijf zn op twee manieren 
EGO rf Ot rk dee ee rx if Og ore r oknageXneg) + 
FOX) ee ee erg) rx df Og sees ox) 

en 
EGOEf Gt rx dees er Eg oee org Xn) 
FOX) ere) Oerd Og oe ee ok) 5 


de eerste manier beschrijft dat men na interpolatie op x ” 


02 n=? 


achtereenvolgens nog Ln en x, als steunpunten erbij neemt, 


mi 





terwijl men in het tweede geval juist eerst Xx, en dan x,_, toer 
voegt. 


Vergelijking van deze twee schrijfwijzen leert : an 


El ponen ke nst Esens are 4) 
WES Pagana hk Eten ee 
n XX „ 
n n-i 


u 
TN 


Zo vinden we voor n Flatt) & f(xjd-f (xp) 


RE 
Merk nu op dat de lineaire interpolatieformule (1.1) een speciaal 
geval blijkt te zijn van (4.1). — 
Men noemt de uitdrukking FlxyoeeesXde n-de gedeelde differentie. 
Een belangrijke eigenschap, die onmiddelijk uit de definitie 
volgt, is dat FXgoee ex) een symmetrische funktie is van zijn 
argumenten d.w.z. zijn waarde verandert niet als men de x!'en 
permuteert. Zo is f(XpsXysX2sXzsXu)d = FfCXz Xu Xi sXorX2) « en 
Derhalve kan men naast (4.2) algemeen schrijven 


mn 


ATEN Ees EE vk AO EREN m1 Amt rn ) 
(u,3) FOgseeesX) = ek 
mk 


Dit stelt ons in staat de gedeelde differenties stelselmatig 
te bepalen. We geven dit schematisch weer in het onderstaande ee 
diagram : 


fz À 
EC Xo aki ska) 


‚) f(x, ) ÉCxXj »X2 1X3) f(xg „XI sX2 2X3 sd 


FCH oN 5% 5 Ad 


rs PERSEE, 

Ne E 
Se 

Pl 


FlHs ss an À 


waarbij men de gedeelde differenties telkens bepaalt uit hun 


onmiddellijke buren linksboven en linksonder. Om de orde te ver- 





UB, 


UC, 


hogen moet men een vanuit f(x,) schuin omhoog lopende rij toe= 
voegen. 

Doordat in het voorgaande niet is geëist dat de steunpunten ge- 
nummerd zijn in de volgorde van grootte, is men bij de praktische 
uitvoering nog vrij wat betreft volgorde awk men de steunpunten 
toevoegt. Ter illustratie geven we een voorbeeld. 


Stel dat we van de funktie f een equidistante tabel hebben, 


zeg voor argumenten 00450 sensi, Als men f(x) wenst te 
kennen voor Xx = 3.274 dan kan men kiezen Rs 3.2,X, = 3.3;X, z Zell, 
x. = 3.5 etc. De op deze manier verkregen variant van (L.1) staat 


3 


bekend als Newtons voorwaartse formule 

Men merkt echter dat naarmate men meer steunpunten neemt, het 
punt x steeds verder van het midden van het interpolatieinterval 
komt te liggen. Volgens (3B) is dat minder gewenst. Om x toch 
steeds ongeveer midden in het interpolatieinterval te hebben, 
kan men kiezen Xx, = 3.2,X, = 3.3,X, = 3.1,X3 7 3-HyX, 5 3.0 etc. 


Dat geeft de formule 
(u.5) £GO=£(3. 2) (xr3.2)E(3. 2,3. 3) t(KT3. ID) (Xx 3. 3) F3. 2,33, DD + 


+(x-3.2)(x-3.3)(x-3.1)£(3.2,3.3,3.1,3.U) +... 
Dit is de formule van Gauss. Merk op dat f(3.2,3.3,3.1)=f(3.1,3.2,3.3); 
zodat men de benodigdel gedeelde) differenties kan aflezen uit een 
schema als (4.4) of (4.7) hieronder. 
Het is te verwachten dat men bij een gelijke interpolatiegraad 
met de formule van Gauss op de aangegeven wijze een nauwkeuriger 
antwoord krijgt dan bij Newton (uit fig.3.5 ziet men dat bij 
vijfde graads interpolatie de fout bij Newton 5x zo groot kan 
zijn als bij Gauss.) 
In het equidistante geval kan men zich werk besparen door niet 


een tafel van gedeelde differenties op te stellen maar van (on- 





(U. 7) 


uD. 


gedeelde) differenties : men definieert recursief 


Kixsh = Elk d 
(1.6) je u 


In een schema ziet dit er als volgt uit: 
Elz) 
AX sx) 


f(x, ) AC sXj »X2 ) 


DS 


f(x) 


5 ACX, sX2 sX3 ) 


ACx, „x4) 


TA 


f(x) 


Men neemt steeds de linkeronderbuur minus de linkerbovenbuur. 

De differenties in de (mt1)-ste kolom van dit schema heten 

m-de orde differenties. 

De formule van Gauss krijgt nu voor het in UB genoemde voorbeeld 


(x = 3.274) de gedaante 
£Gx) = £(3,2) + BAC3.2,3.3) + PPD (3.1,3.2,3. 3) 4 


PrOPPD) (3,1,3.2,3.3,3) + oe 


waarin p = 0.074 (ga na dat dit inderdaad rekenwerk bespaart). 
We bekijken nu hoe men tijdens het verhogen van de orde van de 
interpolatie een indruk kan krijgen van de fout. 

Laat x * x,sXjs--«-. het punt zijn waarin we de funktiewaarde 
moeten benaderen. Veronderstel dat we door interpolatie op 

3 


«X de benadering fr-4 0) hebben gekregen. In plaats van 


02 


n= 


Xn toe te voegen nemen we, in gedachten, x als steunpunt erbij 


m.a.w. we bepalen een n-de graads polynoom Ee zo dat En (xj) = fx) 


voor iz0,...sn-1 en fr) = f(x). 


Er geldt dus (vgl.(H.1)) 


Ex) = E49 + (Xexgdeer ese Cat Elgers get) 


Dit vergelijkend met (3.4) zien we 


se) 8 
Ep see sok) ET voor zekere 5. 
Evenzo geldt 
Mn) 
FCX sere sXpe4? En) Et voor zekere n. 


Als „m) slechts weinig varieert, volgt hieruit dat 
Eg oere sXng®) a EX oere sXnan? en als schatting voor de 


fout gebruiken we 
(4.8) Gek) Geej dees Geer Ê Oger voren) 


Merk op dat (4.8) precies de term is die men bij fj optelt 
om Ela) te krijgen. We stoppen als(4.8) klein genoeg is. 
De aanname dat sm) niet noemenswaardig varieert is essentieel. 


Door te kijken of FCXj oere Xian? voor i = 0,1,2,... slechts weinig 


varieert kan men zich een zeker vertrouwen hieromtrent verschaffen; 


je En) k het bij 
tn Sire voor zekere &j en het zou Di 


(n) 8 nú ze 
toevallig zijn als Eg oeeeoXjan? slechts 


immers EOxgoree sx 
snel variërende f 
langzaam verandert. 

Eigenlijk is dit hetzelfde soort vertrouwen als we in het verhaal 
in 3C moeten hebben. Men kan daarbij weleens bedrogen uitkomen. 
Zo zijn voor f(x) = sin x en Xx, km alle gedeelde differenties 
nul en het proces van ordeverhoging stopt na berekening van f,(x) 


en levert O0 als benadering. 


81. 


82, 


Hoofdstuk 2. Benaderde integratie. 


Inleiding 
b 

Bij de berekening van integralen Sf(t)dt zullen we vaak op 
a 


funkties stuiten waarvan moeilijk of in het geheel niet ex- 
pliciet een primitieve aan te geven is. Om dan toch een benadering 
voor de integraal te vinden kunnen we bijv. de integrand 

door een interbolartepslynsom vervangen en dit polynoom vervolgens 
integreren . Hoe goed zo'n benadering is, zullen we uitgebreid 
onderzoeken. Tevens zal blijken dat we hiervoor het interpolatie- 


polynoom niet expliciet hoeven op te stellen. 


Kwadratuurformules 
We zullen steeds aannemen dat f voldoende vaak continu diffe- 


rentieerbaar is op [a,b]. 


Laten Xp see sXn Ela,bl gegeven zijn. Zij p het Lagrange inter- 
polatiepolynoom van f op de punten Xpo sXje Schrijf f = p+r 
dan is 

b b b 

Sf(x)dx = S[plxldx + fr(x)dx 

a a a 
We weten 


n 
px) = 2 f(x) L Cac) 
hg ER 
Dus 
Db n b tad 
Sp(xldx = EX f(x) S Li (x)dx 
a k=0 a 


b 
De getallen Wie 5 Ni LP ax zijn constanten, onafhankelijk van 


a 


b 
f. Met R= f rlx)dx is dus 
a 


Db n 
Sf(x)dx = À w fx) +R 
a k=0 
n 
We noemen 2 w‚ f(x) de bij X,5.:«sX_ horende interpolatoire 
gen È Ep ì d 
kwadratuurformule voor ff(x)dx. De Xi heten de steunpunten, de 
a 
Wi de gewichten van de kwadratuurformule. 
Stelling 2.1. Is p een polynoom van graad < n, dan is 
b n 
Íp{xidre E ww mla Jd. 
a k=1 K k 


Bewijs : Uit de eenduidigheid van het interpolatiepolynoom (h.1,(2.1)) 
volgt dat p samenvalt met zijn interpolatiepolynoom op de punten 


X 


oor zodat r(x) = 0. 


We kunnen (2.1) ook gebruiken om de gewichten Wi uit te rekenen. 
Neem nl. maar voor p(x) achtereenvolgens leaner ee Dit levert 
een stelsel van n+1 vergelijkingen in de n+1 onbekenden Wi OP: 

Men kan aantonen dat dit stelsel eenduidig oplosbaar is. In de 
hierna te beschouwen voorbeelden zien we die eenduidigheid vanzelf 
(in verband met optredende afrondfouten is het voor grote waarden 


van n niet aanbevelenswaardig de coëfficiënten inderdaad op deze 


wijze te bepalen). 


Voorbeelden 


1. Neem a = Xx, = 0, b=xy=1, hel. Dan is de kwadratuurformule 


Wof(0) + w‚f(1). We weten dat moet gelden 


A 


1 =f 1dxz= wo + wi 


A O 


oi 
1 


2 ie Xdx = Wg-0 + W‚-1 = W‚ 


Poi 
. 


zodat Wy = Wijs 


2. Neem a = Xx, = =Îs Xj 5 0, b = Xx, = 1, n = 2. Op dezelfde 


wijze als boven vinden we 


83. 


3A 


Wo + Wj + W‚, = 2 _ 
Wy + W‚ = 0 
md en 
Wo + WE 5 
” . _ B! _ k 
Hieruit volgt w, = W‚ = FW) Fe 


Om de gewichten van analoge kwadratuurformules op andere inter- 
vallen te vinden, kan men opnieuw deze berekening uitvoeren. Ee 
Dit is echter niet nodig zoals uit de volgende eigenschap blijkt. 
Stelling 2.2 Als 2 w;f(x;) de ie Xe eerXn behorende inter- 


polatoire kwadratuurformule voor f f(x)dx is, dan is 2 Zw;gly;) 


a 
de bij Yys---sy. behorende interpolatoire kwadratuurformule voor 
0 n 
d-c bce-ad 


S gly)dy met Vi Sn Mt 


Se 


Bewijs. Zij a : [a,b] > [e‚dl de transformatie 


alx) = 2 “+ 2 . Kennelijk is a bijectief en er geldt dat 
Pas aCx;). — 
Zij g het interpolatiepolynoom van g op de punten YorrrsYj 


land 


Zij f = ge a. Dan is f = ge a het interpolatiepolynoom van f 


op de punten Xpose sX,° Nu is en 
dr ze b n 
_ d-e ? = d-c d-e 5 
Î sld en zen ze 
Ô s(y)dy ee À f(x)dx ep 0 w‚flx;) en 3 w‚ely;) _ 


Hiermee zijn we klaar. 


Deze eigenschap ziet er ingewikkeld uit, maar er staat in feite niets 
anders dan dat de gewichten evenredig zijn met de lengte van het 
interval en dat de steunpunten in de beide intervallen gelijk- — 


vormig liggen. 


Speciale gevallen ; schatting van de fout. 


Door nulde orde interpolatie van f in het punt m = J(atb) op het 
b 

segment [a,‚bl vinden we de kwadratuurformule voor f f(x)dx 
a 


£ 341 (b-a) f(m) hd 


sB 





Dit noemen we de midpuntregel. In een plaatje ziet het er als 
welet wit 


b 
We benaderen J f(x)dx door het ge- 


Ee a 
Dx EN arceerde oppervlak. Om een schatting van 


de restterm te geven ontwikkelen we f rond 


m in een Taylorreeks 


NS 
Ss (ERI Bloed = EC) + Gem) fm) + EE EERE) 
a m b 
met E(x) tussen x en m. We integreren nu (3.2) 
Er komt 

b b leni? 

S fax = (b-adflm) + Sf EE FME dx 

a a 


(alhoewel we niet weten hoe E(x) van x afhangt blijkt uit (3.2) 


onmiddellijk dat FY(E(X)) een continue functie van Xx is; 


(x-m) ? 
2 
dus integreerbaar; overigens is het ook niet moeilijk aan te 
tonen dat f"(E(x)) een continue functie van x is als we nog 

definieèren Em) = m, maar dat hebben we hier niet nodig). 
Wegens (x-m)? > 0 ligt de integrand in het rechterlid tussen 


(x-m)?2 min f" en (x-m)? max £'‚, dus is er een E E(a,b) zodat de 


integraal rechts gelijk is aan 


b 2 _a3 
PUEEN F emd ax - Bal grs) asE<b 
a 
Zodat 
Fi (b-a)® 
(3.3) S flx)dx = (b-a)flm) + == fY(E) a<EsS<b 


a 
Een tweede geval dat we beschouwen is dat f benaderd wordt door 


a, Xx, = b. Dit hebben 


lineaire interpolatie op de punten x, 


we in voorbeeld 1 van 82 gedaan met a = 0 b = 1. We hebben hier 
ook de gewichten uitgerekend : Wo = Wj = à° 
b 
Met (2.2) vinden we algemeen de kwadratuurformule van Jf(x)dx 
a 
(nat DE drh: de ED 
5 D 0 1 


Cen 


Dat deze formule de trapeziumregel 
heet, is uit nevenstaande figuur 


D 


direkt duidelijk. We benaderen ff(x)dx 
a 


weer door het gearceerde oppervlak 
a b 


Om de restterm te bepalen schrijven we 
u 5 
roo = EEEN Gek aex) 


Omdat r(x) tekenvast is, vindt men geheel analoog aan het boven- 


staande 


â 
(3.5) Re {be) £"(E) at LE 


Tenslotte behandelen we de regel van Simpson. Hierbij is 

Xo = a, xj = Flatb)sx, = b 3 we benaderen f met het tweede orde 
interpolatiepolynoom op X,» Xj» X,- Op analoge wijze als bij de 
trapeziumregel vinden we de gewichten, nu met voorbeeld 2 van 
52 waar we het speciale geval a ="1, b = 1 bekeken hebben ; met 


b 


(2.2) vinden we de kwadratuurformule voor Sf f(x)dx 
a 


(Je B 


wi 


[ECD + FC) + FC) 


waarbij h = Ee 


2 
Bij Simpson (waar men op grond van stelling (2.1) slechts exacte 
uitkomsten voor hoogstens tweede graads polynomen mag verwachten), 
constateert men wellicht enigszins verrast, dat men op het segment 


[-1,1] ook nog de exacte waarde van de integraal voor derde graads 


polynomen krijgt. Immers 
k 1 

f x'axe 0 = 2l-1)* + 

-Ì 3 


k o3 1 4% 
3 0 + 3 Ls 


Dus integreert Simpson 3e graads polynomen exact op elk interval 


omdat de in het bewijs van stelling (2.2) gebruikte substitutie 


derde graads polynomen in derde graads polynomen overvoert 


b ttt 
Er geldt R = jE gE (xrx,d(xrx, )(x-x,)dx. Helaas is 
Ee : 


(xexy) (xrx,)(xex,) niet tekenvast en deze misère komt men vaak 
tegen als men kwadratuur resttermen wil bepalen. We ondervangen 
dit met een truc. Zij q het derde orde interpolatiepolynoom van f 


zodat q(x;) = flx;), Ì = 0,1,2 en q'xj) = f'Cx,)s Dan geldt 


(4) 
f(x) = QX) + EE (xk, ICRm 1E Cxex, ) 
Dus 
b b b(H) 
fElxhax = Jqlxjde + je E) (HR) kek) Cres IdX 
a a a 


omdat q een derde graads polynoom is, geldt 
D 
Keken = EZ wjalx;) = 2 w‚flx‚) 
volgens de definitie van q. Bijgevolg 
D b_(4) 
JE Khel =Ew;flx) + ek (x=) xk dP (Xexy dx. 
a 


Hierin is Cot Wed Cor) nu wel tekenvast. Derhalve geldt 


CH) 


R = Ere) , En Tam zj) Cun, dax 
a 
oftewel 
En 5 
(87) R=- dte) f®cE) ae 
90° 2 
Opmerking. Als een kwadratuurformule voor een interval [a,b] een 


restterm van de gedaante pf le) heeft, dan heeft de volgens 


stelling (2.2) bijbehorende formule voor het interval (c‚d) de 
ee 


restterm Es) (E). Dit blijkt door toepassen van de sub- 


b-a 
stitutie gebruikt in het bewijs van stelling (2.2). Vandaar dat 


men in resttermen steeds tegenkomt neleke, 


SU, 


Gepenekeende kwadratuurformules. 


Een kwadratuurformule is zelden nauwkeurig genoeg om op een 
interval van enige lengte een redelijke benadering van de 
integraal op te leveren. Men verdeelt dan het interval in een 
aantal delen en En die kwadratuurformule toe op elk deel- 
interval. Een op deze wijze gevormde kwadratuurformule noemt 
men een gerepeteerde kwadratuurformule. 

Verdeel het interval [a;,bl in n stukken, ter grootte h. Noem 


+ ih. De trapeziumregel op een deelinterval luidt 


X 
X 
es A 0: EN 3 en 
E Ê(x)dx = zl fx) + fx)! 15 hef (5) 


2-4 


1 
b dre 
LER = hi BEKO) +ECK DE tf Ox FEE OD Eph Rr 


omdat f' een continue functie is is het gemiddelde van de 
waarden ENE) sere Ê TCE) weer een waarde van f" 3; er is dus 


een & Ela,b) zodat 


n 
B fUlE) == nfE), 


1=1 
Met nh = b-a vinden we zo 
b 
(u. 1) hae = hf )tfO tf Dt. tE Or tEE 
„-Â (b-a)n?f"(E) 
12 


De aldus verkregen kwadratuurformule heet de gerepeteerde 
trapeziumregel. 

Ook de formule van Simpson kunnen we herhaald toepassen. We 
doen dit op de intervallen [XysXpiagl is hier dus even). We 


vinden zo de gerepeteerde formule van Simpson. 








: | 
(4.2) JEGOAx = BEG )EBE(K DECK DEED one + BEG D+EO)I 
a 
was (saint) 
ig braomn'el Pe) 


Ga dit zelf na. Let er bij de afleiding van de restterm op dat 


het aantal intervallen waarin (a,b) verdeeld is, in is. 


Opmerkingen. 1. Men ziet uit (h.1) en (h,2)-dat als f!' resp. 
£() weinig varieert op het integratieinterval in kwestie, 
halvering van h bij de trapeziumregel een tx en bij Simpson 
een 16x zo goed resultaat geeft. Dit is niet in strijd met de 
indruk die men uit (3.5) en (3.7) zou kunnen krijgen, nl. aes 
dat dit 8x en 32x is: immers, halveert men het integratie inter- 
val dan moet men, om de oorspronkelijke integraal te berekenen 
de kwadratuurformule ook nog op de andere helft toepassen. 

2. Het is duidelijk dat voor h > 0 het resultaat 
van de beide kwadratuurformules nadert tot de ware integraal. 
In plaats hiervan zou men ook eens kunnen proberen door een- 
il ige toepassing van steeds hoger orde interpolatoire kwadratuur- 
formules steeds nauwkeuriger resultaten te bereiken, 
Hier geldt echter dat convergentie precies dan plaats vindt 
als de sommen der absolute waarden der gewichten per zo ver- 
kregen kwadratuurformule een gemeenschappelijke bovengrens 
hebben. Dit nu blijkt niet het geval te zijn in het voor de 
hand liggende geval dat men de steunpunten equidistant kiest 
(de zgn. Newton-Cotes formules). Hiervan worden de gewichten 
op den duur sterk positief en negatief. Het gaat echter wel 
geed bij de formules van Gauss, die we in de volgende para- 


graaf summier zullen behandelen. 





ge: 


Voorbeelden. In onderstaand tabelletje geven we de restterm weer 


Gauss formules. 


De algemene gedaante van een kwadratuurformule luidt 


b N 
(5.1) ÍElxddr == -Z w; f(x) Eeen 
a 1=0 


In het voorgaande hebben we de steunpunten x, steeds als vaste 
gegevens beschouwd. Dit hoeft echter niet. Dan blijkt dat (5.1) 
2N+2 vrij te kiezen parameters heeft, en het ligt voor de hand has 
deze zo zimstie mogelijk te kiezen. We vermoeden dat het gunstig 

is als polynomen van zo hoog mogelijke graad exact worden ge- 
integreerd, en men kan bewijzen dat dit inderdaad zo is. Door 

te eisen dat A exact worden geintegreerd, ontstaat 

een stelsel van 2N+2 niet-lineaire vergelijkingen in de 2N+2 
onbekenden w;‚x,- Men kan aantonen dat dit stelsel eenduidig 

oplosbaar is, zie b.v. de syllabus Numerieke Analyse 1. 

Steunpunten en gewichten kan men in tabellen-boeken vinden. Bijv. 

(zeer aanbevelenswaardig) Abramowitz-Stegun: Handbook of Mathematical 


Functions. Dover publ. 1965 ($4 voor 1043 pag.) pag. 916. We — 


vermelden slechts de restterm : er geldt voor n+1 punts Gauss: 


2n+3 n+3 


_ (bead Pt Sn+1!* C2n+2) otsibeat? C2n+2) 
KE (2n+3)((2n+2)!)3 d (5) u2nt2ron+2)! 5 (6) 


(dit laatste door Stirling toe te passen), a < 5E Sb. 


Ti 
als men ter berekening van fsin(x)dx een zeker aantal punten inves- 
0 n 
teert in trap., Simpson en 5-punts Gauss. 
trap Simpson 5-punts Gauss 
5 punts „iej ged =d 8 Leben? \ ik 
9 punts =3 5402 A mn 1.2,,710 


17 punts „Dg 8 2309 del, 48 


86, 


Opvallend is dat 5 punts Gauss enkelvoudig toegepast reeds zo'n 
goed resultaat geeft, alsmede dat halvering van het interval bij 


9 punts Gauss de nauwkeurigheid zo snel doet toenemen. 


Glijdende stap 


Als we een functie f over een interval [a,b] numeriek willen 
integreren met een gegeven tolerantiezep e‚ dan kunnen we b.v, 

de gerepeteerde regel van Simpson gebruiken. Hierbij vallen 

direct twee bezwaren op. 


Ch) 


1. Om de stapgrootte te bepalen, moet men iets weten over f e 


(u) 


Schattingen van f zijn echter veelal moeilijk te verkrijgen. 
Ze Als f een verloop als in de figuur 
heeft, dan zal de stapgrootte be- 
paald worden door de steilheid in de 
pieken. Deze stapgrootte is dan op 


de rest van het interval veel te 


fijn. Er wordt hier dus veel over- 





bodig werk gedaan. 
b 
In het volgende zullen we een algorithme bepalen, dat Jf(x)dx met 
. a 
zekere nauwkeurigheid oplevert, zonder dat hier schattingen van 


£€) voor nodig zijn, en zodat het bezwaar onder 2 wordt onder- 
vangen. Men zij echter reeds nu gewaarschuwd dat de methode niet 
geheel waterdicht is; zie de opmerking aan het eind van deze 
paragraaf. 

Laat I,(a,h) de regel van Simpson zijn toegepast op het interval 
[a,a+2hl , en zij 


dij 


Ij(a,h) = Io(asd) + Ig(ath,g 


de eenmaal gerepeteerde regel van Simpson. Dan is volgens (3.7) 


a+2h 


(6.1) 1 £Gdx= Ilam) - do hele) 
a 
Evenzo is volgens (4.2) 
at2h 
4 
(642) Of EGOaxE Ia ie do nief 


Als we nu eens aannemen dat £() op het beschouwde interval 


(u) 


weinig varieert, dan kunnen we de term f£ uit (6.1) en (6.2) 


elimineren. We vinden zo 


at ?2h 1 
(6.3) S f(x)dx - I,(a,h) = 15 CT, LAAN) — Tod 
a 


We hebben dus een mogelijkheid gevonden om de fout, die I, 
laat, te schatten (de waarden van I, en 1, kent men! ). 
We komen nu tot de volgende werkwijze : Begin met h = 3(b-a) 


en bepaal TI, en E Kijk met (6.3) of de fout al kleiner dan 


1 . 


at2h 
e is. Zo ja, dan accepteren we I, als benadering van S ECX)dX = 
b a 


z= Sf(x)dx. Zo nee, dan halveren we h, en we herhalen de procedure; 
a 

maar nu met tolerantie e/2,voor de linker helft van [a,bl. 

Dit proces herhalen we zo lang totdat, zeg na n keer, voldaan 


is aan de eis dat 


(Tj (ash) = Ilan) < se 
at2h 

Dan accepteren we I,(a,h) als benadering van { f(x)dx. 
Nu we la,a,) afgewerkt hebben, gaan we à 
verder met [a sa,l- Dit is de rechter- 
helft die over was van de laatste halvering, nl. die van 
[a,a,l. Deze helft werken we met de glijdende stap methode af. 
We hebben dan de integraal op [a,a,l benaderd. Merk nu op dat 


[a‚a,l ontstaan is door halvering van [a,a,l. Als boven gaan 


we nu de rechterhelft doen, dit is [a,sa,l, enzovoort. 


Opmerking. De glijdende stap algorithme is niet waterdicht, 
Cu) 


tenzij f over het hele interval weinig varieert. Dit laatste 


is bij ingewikkelde functies natuurlijk niet het geval. En’ 
zodra g6 nulpunten heeft kan het natuurlijk helemaal misgaan 
Toch valt het in de praktijk wel mee. Stel nl. eens dat op het 


hie Alt. 


aan de beurt zijnde deelintervalletje zo'n nulpunt van f 
De foutschattingsformule (6.3) levert dan bijv. een 10x te 
kleine waarde op, zodat op dit deelinterval de beschikbare 
tolerantie ver overschreden wordt. Deze beschikbare tolerantie 
is echter doorgaans een klein deel van de totale e,‚ zodat deze 
overschrijding op het geheel weinig van invloed is. 

Maar het kan echt goed misgaan 
Nevenstaande prent stelt de 


grafiek van een functie voor 


waarvan de integraal # 0 is, 





terwijl T,b=l,4 s I,l-d,4) = 0, 
zodat het proces na de eerste slag al stopt en O0 als uitkomst 


geeft. 


8. 


Hoofdstuk 3 Numerieke behandeling van differentiaalvergelijkingen. 


Oplosmethoden. 


Uitgaande van numerieke integratie heeft men een aantal 


methodes ontwikkeld om de oplossing van een beginwaardeprobleem 


Sateen 


op een segment [t‚‚t,‚+T] te benaderen. Enkele van deze methodes 
zullen we in dit hoofdstuk behandelen. Existentie en uniciteit 
van de (exakte) oplossing staan hier niet ter discussie; we 
zullen steeds veronderstellen dat er precies één oplossing is. 
Een benadering van de oplossing bestaat uit een "tabel! met 
benaderde waarden van de oplossing in een aantal punten Ee van 
het segment. Aangezien t vaak als tijd geïnterpreteerd wordt, 
zullen we meestal spreken over de tijdstippen t: Deze tijd- 
stippen kiezen we equidistant, dus t, 5 to + nh ; h heet de 
stapgrootte. 


De te behandelen oplosmethoden berusten op de opmerking dat 


t 
n 
(1.2) dB atf Plts Gadet 
E 
n-1 


waarbij Xx’ de exacte oplossing van (1.1) is en Xn de waarde 
hiervan op het tijdstip tn: : 


Nu geldt voor een willekeurige differentieerbare functie g bijv. 


ie 
n 2 
2 h° 
bla) Ki gitjdt = Hett l + 5 BE dst SE S B 
n-1 


(integratie van het O-de orde interpolatiepolynoom in Brei” 


Door toepassing op (1.2) komt er 


2 tt 
(1.4) É En PRE ot d- Xx Ede t SE Et 





Analoog krijgt men een kwadratuurformule als erachter aangeduid 


_ * h° * e . . 
EED, KS Beg td Ee gd WE) KU ge. Anterp. pol. int) 
(1.6) EB ú + in [fl + fl] - hi ede ) (trapeziumregel) 
fi n n=1 « n n= 1 12 n 8 
Hiesta is mts sh. 
n n° Sm 


Door in deze formules de termen met 5 weg te laten krijgt men 


relaties van de gedaante 


Ci. 7 KH + hf (Euler (forward)) 

n n=-1 n-1 
(1.8) Ee + hf (Euler backward) 

n n=1 

ek 1 . 

(18) en Mon 3h[f, + Eed (trapezium) 
waarin in = ECE Xx): 
Men ziet dat deze relaties (1.7) - (1.9) uitgaande van een getal 
Ky Sen wij Ky sXgrXgoree definiëren, en men hoopt dat Xn onder 


omstandigheden een redelijke benadering voor Xn is (dergelijke 

betrekkingen noemt men recurrente betrekkingen). Deze methoden 

om de differentiaalvergelijking benaderd op te lossen duiden we 
aan met de erachter staande namen. 


2 Li 
Omdat het net is alsof we in (1.7) t/m (1.9) de termen a xXx (E) 


n 
uit (1.4) en analoog die uit (1.5) en (1.6) verwaarloosd hebben, 
noemen we deze termen de locale (= plaatselijke) discretisatie 
fout. In de literatuur wordt deze benaming ook vaak gegeven aan 
bedoelde termen gedeeld door n; in verband hiermee noemt men de 
methoden (1.7) t/m (1,9) van de orde wesp. 1, 1 en 2. 

We merken nog op dat de methode van Euler er grafisch blijkbaar 
op neer komt dat men in elk punt CE’ een stapje doet in 


de richting van het richtingsveld ter plaatse, terwijl men bij 


Euler backward uitgaande van CE) een punt (tsx) zoekt 


waar het richtingslijntje door (t > Poikt, 


n-1’*n-1 
In £1.7) komt x, alleen in het linkerlid voor en is dus ge- 


makkelijk te berekenen, als men x al heeft; (1.7) is dan ook 


n=1 
een voorbeeld van een expliciete methode. De methodes (1.8) en 
(1.9) zijn daarentegen impliciet: ze geven een vergelijking waar- 
Hit Xx, moet worden opgelost. Veronderstel dat de vergelijking 


oplosbaar is, dan kan men als volgt te werk gaan om deze oplossing 


te benaderen (zie ook hoofdstuk 4). 


De recurrente betrekkingen (1.7) - (1.9) zijn allemaal van de 
vorm 
CT 0) Br Hat h(a £CE, sx) + B FCE gep) 


Kortweg 


X 
n 


dx) 
(0) 


Schat nu Xx, en geef deze schatting aan met Xx, 3 men kan bijv. 


ak = Xp-j remen. Definieer vervolgens 
(IJ Cad Ke 
Xn = dx ) voor Ì 21. 
Dan geldt 
CER — ES 
(1.11) xj k ® dx, ) ò Cx) 
Ci) 
_ 1 pn 
z= Ò'(E) xn xj) 
voor zekere E tussen gr en Xx 
Als Ens ‚ de afgeleide van f naar de tweede variabele, begrensd 


is, dan geldt voor voldoend kleine h 


[é' CE) = [ha EJA < constante < 1. 


Ci), 


En m.b.v. (1.11) zien we dat voor voldoend kleine h de rij (x 
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naar X, convergeert. 


Foutengedrag. 


° Ee * 
We bestuderen nu de verschilrij le}; e= Xj 7 Xj» van een 


benaderde en een ware oplossing van een diff. verg. De algemene 
theorie hiervoor is niet eenvoudig. Wat er kan gebeuren wordt 


echter al heel aardig geïllustreerd door de simpele dvgl. 


Cet x'(t) = pxlt) + q(t), x(0) = x, 


voor t € [O,T]. We bekijken eens de methode (1.7), met t‚ = nh. 
(2.2). X = X + hf 
Hiervan (1.4) aftrekken geeft 
n? , 
(2.3) Bis Bt CL + hp} Bd? (NB: Es se ger (GE dls 


Schrijf 1 + hp = 0. Als we de recurrente betrekking in (2.3) uit- 


schrijven krijgen we 


(2.4) es Sk GE HE H 

= 2 n-1 

= Ö, + oÔ + 0 Ön-2 + + 0 ô, + 0e, 
Aangezien x, gegeven is geldt e, = 0. 


We zijn er nu uiteraard in geïnteresseerd wat er gebeurt voor 
kleine waarden van h. Maar men moet wel oppassen, want als we 

h zonder meer naar 0 laten gaan (d.w.z. n vasthouden) dan gaat 
het punt t‚ naar O0 en geeft e, de fout aan in een punt dat naar 

0 gaat en dat is niet zo interessant. We willen graag weten wat 
de fout doet in een vast punt t als we h naar 0 laten gaan, d.w.z 


we zullen moeten toestaan dat n stijgt als h daalt. Het ligt 


daarom ook min of meer voor de hand dat we voor kleine h het — 
rechterlid van (2.4) proberen te vervangen door een integraal. 

ki 2 - 
ein(1+ph) E „Ph+0(h ) 


Daartoe schrijven we 0 = zodat wegens 


hm S$ T 
(2.5) es ePhm + O(h) ; 


Gebruiken we verder dat de lokale discretisatiefout van Euler 


2 t 
8, ze EE TED is, dan vinden we voor (2.4), weer wegens 
hn < T 
2 1 "u tr N 
e= Br ce) + PP En) + PPP ED tel + OCH) 
(2.6) pn 
2 p(t st ) u p(t RE ” ) u 
== le a CE dte B MEE ef LN Rn 
DOE =d 
ek ie ME dbs el + O(h®) 


Nu geldt voor een willekeurige differentieerbare funktie f dat 


Es t 
Eh IN ie 
hed ed fls)ds + O(h?) (ga na), zodat hE f(t‚) = f f(s)ds + O(h) 
Es 0 0 In 
dsl 


Door dit toe te passen op (2.6) vinden we dat de totale fout in 
de benaderde oplossing op het tijdstip t gelijk is aan 

nt p(t-s)e" é 
WEP = 7 Ä e x (s)ds + 0 (h*) 
Deze formule laat allereerst zien dat de benaderde oplossingen van 
(2.1) uniform naar de exakte oplossing convergeren als h naar nul 
gaat. 


Belangrijker echter voor de toepassingen is, dat de formule be- 


wijst dat bij Euler de fout, op hogere orde termen na, evenredig 


is met de stapgrootte h; vandaar de benaming eerste orde proces 
in 61. Als men dus de waarde van de oplossing voor een zekere t 
eerst eens benadert met behulp van een zekere h en daarna nog 
eens met de halve h, dan zal de fout in deze laatste uitkomst 
ongeveer gelijk zijn aan het verschil van de twee uitkomsten. 
Hiermee heeft men dan een schatting voor de bereikte nauwkeurig- 
heid. Blijkt nu dat de fout (bij Jh) tienmaal zo groot is als men 
wenst, dan is het duidelijk dat men het nog maar eens moet pro- 


beren met stapgrootte De 


Het is nu ook duidelijk wat men bij een tweede orde proces, zoals 
(2.8) mag verwachten : de totale fout op een tijdstip t is even- 
redig met h? en de fout bij stapgrootte Jh is driemaal zo klein 


als het verschil tussen de uitkomsten bij h en 5h.- 


Men moet een beetje voorzichtig zijn met dit soort uitspraken. 

Als gn niet tekenvast is, kan het gebeuren dat voor een zekere 

t de integraal in (2.7) nul is. De term O(h?) bepaalt in dat 

geval de fout en de fout zal dus wel niet gehalveerd worden, 

als men de stapgrootte halveert. De fout in zo'n punt is dan 

juist erg klein, kleiner dan in de omringende punten. Voor verre- 
weg de meeste punten t zal de integraal echter niet nul zijn, zodat 
men, als men de oplossingen tabuleert en twee zulke tabellen; 


-— de een verkregen met h, de andere met jh - naast elkaar legt; 


toch een goede indruk krijgt van de fout. 
Hoe dit in de praktijk uitvalt toont het volgende voorbeeld. 


Voorbeeld. 
Onderstaande tabel geeft voor de aangegeven waarden van t de ware 
oplossing en de fout in benaderde oplossingen,verkregen op de 


aangegeven wijze,van de dvgl. 


dx 
dt 


ols 


= Xt COSCE), Ky == 
met als oplossing 


teh = Â Gine) = BEELEN. 


De homogene dvgl. z = Xx heeft als oplossingen x(t) = C en 


zodat met name fouten aan het begin gemaakt, nogal aangroeien. 


Men ziet de fout bij Euler forward heel aardig evenredig is met 


h, (behalve daar waar de fout zowat 0 is) en dat Euler backward 


83. 


geen kwaliteitsverschil geeft met Euler forwards. 


+00. 3 
+00.6 
+00,9 
+01.2 
+01.5 
+01. 
+02, 
+02. 
+02. 
+03. 
+03. 
+03. 
+03, 
+04, 
+OU, 
+OU, 
+05. 
+05. 
+05, 
+06. 


OENE WO AIEE 


x(t) 


=l 


-0. 


+0 


+0. 
+0, 
+0. 


+0 


+0. 
+0. 
+0. 


+0 


+0. 
+0, 
0. 
=). 


-Û 
= 0 
=d 


„0, 


ade, 


‚330 
130 
‚081 
285 
463 
601 
‚684 
706 
666 
566 
‚415 
227 
019 
191 
383 
‚542 
„652 
„704 
693 
„620 


fout voor Euler forward 
SO 


h = 0.01 


-000.0007 
„000. 0014 
-000.0019 
-000.0023 
-000.0025 
-000.0024 
-000.0021 
-=000.0017 
-000.0010 
-000.0003 
+000.0005 
+000.0012 
+000.0018 
+000.0023 
+000.0026 
+000.0027 
+000.0026 
+000.0024 
+000.0020 
+000.0015 


Foutvoortplanting. 


h = 0.02 


-000.0015 
-000.0028 
-000.0039 
-000.0046 
-000.0049 
-000.0048 
-000.0042 
-000.0033 
-000.0020 
-000.0005 
+000.0010 
+000.0025 
+000.0038 
+000. 0049 
+000.0056 
+000.0059 
+000.0059 
+000. 0056 
+000.0051 
+000. 0045 


h = 0.10 


000071 
-00.0136 
-00.0188 
00222 
0D, BZS 
„-00,0226 
-00.,0195 
„00.044 
„000076 
+00.0002 
+00.0087 
+00,0171 
+00.0251 
+00,0322 
+00.0384 
+00.0435 
+00.0481 
+00.0528 
+00.0587 
+00.0670 


backward 


A 


h = 0.10 


+00,0077 
+00.0147 
+00.0205 
+00.0246 
+00,0266 
+00.,0265 
+00,0243 
+00.0202 
+00.0149 
+00.0090 
+00.0033 
-00.0014 
„00. 0041 
-00.0039 
+00.0001 
+00.0089 
+00,0235 
+00, 0449 
+00.0747 
+00,1147 


We bekijken nog eens de differentiaalvergelijking uit 82. 


Formule (2.4) laat zien dat bij de k° stap gemaakte fout Ö 


zieh bij de n 


n= 
stap als een go 


ae maal zo grote fout manifesteert. 


Of ook: de op het tijdstip t gemaakte fout ê_ blijkt op het 


p(t-t) bh 


tijdstip t al e + O(h) maal zo groot te zijn geworden (vgl.2.5) 


Ter vergelijking belasten we eens de oplossing van de differen- 
tiaalvergelijking ten tijde t met een fout ê. We bekijken dan 
de oplossing X met x(t) = x°(Ì) + 8. Nu heeft het verschil van 


twee oplossingen van (2.1) de gedaante C.ePt zodat 


kt) -_ x(t) == geP(E-t), 


p(t-t) 


Deze fout is dus met een factor e aangegroeid. Dit resul- 


taat lijkt verrassend veel op het hierboven verkregene. 


Globaal sprekend kunnen we dus zeggen dat bij (2.2) het maken 
van een (discretisatie- of afrond-) fout ê bij t ongeveer neer- 
komt op het overstappen op de oplossing x(t) als boven, of ook 
dat fouten zich ongeveer voortplanten als oplossing van fe: 
homogene diff.vgl. Men kan aantonen dat dit gedrag bij andere 
(nette) methoden ook opgaat. Dit geldt echter alleen maar onder 
de aanname dat h klein is, en de vraag wanneer h klein is kan 
hier onaangenaam uitpakken. 


We bekijken de differentiaalvergelijking 





be 1) le „xtcos(t/20) 
De oplossing hiervan met x(0) = x, is 

_ _ 400, _-t , KOO EE ke Ee 
Ca 2À REEN B ER, Tore + Toi ces zo * zo Sin 50) 


Men ziet dat de eerste term in het rechterlid nogal snel uit- 


400 


05 Toi en dat de oplossing zich daar- 


sterft (of 0 is, als nl. x 
6 . - 4 t 1 ei t 
na gedraayt als de langzaam golvende functie cos 50 + Dg SIN 57 


die telkens als t met © 60 toeneemt, een halve periode doorloopt. 


È 


Derhalve, als men deze functie bemonstert voor t = 0,3,6,9,... 
krijgt men een uitstekend beeld van de functie; we kunnen ook 
zeggen dat deze funktie zich met stapgrootte 3 uitstekend laat be- 
tE 


schrijven; in tegenstelling tot funkties als sin(t) of e Tt 


verband hiermee lijkt bijv. h = 3 dan ook een heel redelijke 
stapgrootte voor numerieke oplossing van de dvgl. En inderdaad, 


de locale discretisatiefout bij (1.7) is nu 


LOO _-t _ 400 
"Tore 


ie 


EN 


* 20 20 


ê keye Ae 


Ca. 8) 5 


(cos 5 


oo 


En 
401’ 00 


sol 


en dit is hoogstens © 1% van de oplossing als de kopterm een- 
maal voldoende uitgestorven is. 

De locale discretisatiefout is dus redelijk klein in vergelijking 
met de waarden van de oplossing; bovendien sterven de oplossingen 
van de homogene dvgl. x!' = -x snel uit, zodat we verwachten dat 
ook de effecten der locale discretisatiefouten (die zich immers 
voor kleine h ongeveer als oplossingen van de homogene vgl. ge- 
dragen)snel uitsterven, waardoor dus de totale fout bij een 

niet te kleine t voor het grootste deel zou bestaan uit de 
effecten van de locale discretisatiefouten in de laatste paar 
tijdstappen ende discretisatiefouten uit het begin nauwe- 

lijks meer bijdragen. 


Bekijken we nu echter (2.4). Hierin is o = 1-h = -2, zodat 


(3,4) ë sé, — A + Uê = 80 t zen È 2 , 


en men ziet dat juist het effect van ö, en de overige "vroege! 
8's buitensporig is. 
Hier heeft men dus een typisch geval van slechte foutvoortplan- 


ting, van kleine oorzaken die grote gevolgen hebben. In zo'n 


geval spreekt de numericus van instabiliteit, en hij noemt een 


a 








proces stabiel wanneer kleine fouten onderweg ook maar tot 

kleine fouten in de uitkomst leiden (interpretatie van 'klein' 
van geval tot geval bepalen). 

Blijkbaar was de gekozen waarde van h nog niet klein genoeg om 
ervoor te zorgen dat de fouten zich inderdaad ongeveer voort- 
planten als oplossingen van de homogene dvgl., en inderdaad: het 
is ook wel duidelijk dat h = 3 een niet erg geschikte enkel 
is ten aanzien van de homogene differentiaalvergelijking, omdat 


de funktie 8 niet erg goed beschreven wordt met stapgrootte 3. 


Nu is het niet bij alle methoden zo dat fouten opeens gaan aan- 
groeien wanneer h eigenlijk te groot is om de oplossingen van de 
homogene differentiaalvergelijking te beschrijven. 


Bekijken we (1.8) eens. In plaats van (2.3) vinden we nu 


| n BE en 
{8e B) BS 8 en-1 * HP es bk NEN 
oftewel 
(3.6) es 08 ge met oz a 

n n n=-1 1-hp 
zodat we analoog aan (2.4) krijgen: 
hé 2 n 

(3.7) a, ej Pace t 0 ö, 
Echter is in ons voorbeeld P = -1, zodat er nu met h = 3 komt 
oz z. Weliswaar is weer niet, zoals in (2,8), of = em 
tinie oe Eri en e” tm PN GD, maar het belangrijkste 


is dat fouten nu wel uitsterven en snel ook. Methode (AB) 


laat derhalve wel de gewenste h toe. Analoog (1.9). 


Resumerena kunnen we zeggen dat we in onze ervaringen in deze 
paragraaf tot dusverre de volgende ingrediënten hebben aange. 


troffen 


-— men beschouwt een stapgrootte h als redelijk wanneer de ware 
oplossing van de differentiaalvergelijking behoorlijk beschre- 
ven wordt door zijn waarden in de punten ih, i = 0,1,23...5; 
dit betekent dat de locale discretisatiefout dan klein is ten 
opzichte van de funktiewaarden in hun totaliteit (niet ten op- 
zichte van elke funktiewaarde afzonderlijk, want zo'n waarde 
kan wel 0 zijn). 

=- wanneer de oplossingen van de homogene differentiaalvergelijking 
snel uitsterven ten opzichte van het gedrag van de oplossingen 
van de hele differentiaalvergelijking dan beschrijft boven- 
staande redelijke stapgrootte de oplossingen van de homogene 
differentiaalvergelijking niet goed. 


- dit veroorzaakt bij sommige methoden slechte foutvoortplanting; 


bij andere niet. De laatsten zijn dan verkieslijk. 


Differentiaalvergelijkingen van het beschreven karakter noemt men 


stijf. 


We bekijken ook nog eens grafisch wat er gebeurt. In fig. 3.1 
is het richtingsveld van (3.1) ge- EN 

schetst (NB schaal in horizontale | a En 
richting 20 x zo klein als in ver- 
ticale). De richtingscoêfficiënt in 


(t‚x) is in werkelijkheid -xtcos 55 





maar in onze tekening door de ge-- 
kozen schalen 20 maal zo groot, zo- 
dat reeds op 0.1 boven en onder de grafiek de getekende richtings- 
coëfficiënt + 2 is. Het veld staat dus bijna overal vrijwel 
verticaal, behalve in de onmiddelijke omgeving van de grafiek van 


gt) = cos En ‚ Dus elke oplossing, waar ook gestart, begeeft zich 











onmiddelijk naar g,‚ en blijft daar voortaan dicht bij de buurt. 


In fig. 3.2 ziet men Euler weergegeven. Zodra men door de dis- 
ceretisatiefout iets te ver van de grafiek van g vandaan komt 


valt men in de klauwen van het verticale richtingsveld. In 











fig. 3.2 fig.3.3 


fig. 3.3 staat Euler backward . Gezien de grafische interpretatie 


hiervan kan men niet ver van de grafiek van g weg komen. 


Tenslotte merken we nog op ook bij Euler backward de fout on= 
geveer evenredig is met h, ook al is deze zo groot dat de op- 


lossingen van de homogene differentiaalvergelijking er niet goed 





mee beschreven kunnen worden. Als nl. 0 = Er < 1 dan doen in 
(3.7) slechts o 8,0? 8 rea 8 met een betrekkelijk kleine 
n nel n=k a8 
k mee. Hiervoor is &, zowat constant = Ô ‚ dus e = ET = 
2 Lj ie Ad E S ú 
ERE _h 4 
hp * (5) = ie (E). 


Een en ander wordt fraai vertoond in onderstaand voorbeeld. 
(merk op dat hierin voor Euler backward de evenredigheid van de 
fout met h al mooi opgaat alhoewel o niet eens zo erg klein is 


nl. og = OD, resp. 0.25). 


— U — 
Voorbeeld 
Onderstaande tabel geeft voor de gegeven waarden van t de 
ware oplossing en de fout in benaderde oplossingen verkregen 


op de aangegeven wijze voor de differentiaalvergelijking. 


Eon ee ss 
qe ° Tx + eoslt/20),x, = 0 
met als oplossing 
* _ 400 EE et A 
x(t) = TOT [eos zo * Sn ein 5 el] 


De vergelijking vertoont stijfheidsverschijnselen (snel uit- 
stervende oplossing van de homogene differentiaalvergelijking). 
Euler forwards heeft achtereenvolgens 0 = 5 go = =i enoe= =-2, 
en men ziet dat de grote beginfout uitsterft (h = 1.5), constant 
blijft op het teken na (h = 2), resp. opgeblazen wordt (h = 3; 
-2.625'+05 leze men als -2.625 x 10. 

0.25, en 


Euler backwards heeft achtereenvolgens o = O.b en og = 


gedraagt zich uitstekend. Men ziet dat ook nu de fout weer even- 


redig is met h. 





Euler forward backward 
RN EE EN 
Ë eit Ms dak Hs h = 3 hsl he 8 

+06.0 +0.965 —00.0582 +01.0049 -—-3.999“+00 —-00.0249. -00.0633 
+12.0 +0.851 =—00.0023 —00.9978 =—1.602”+01 =00.0023 —00.0072 
+18.0 +0.659 +00.,.0010 +01.0017 —-6.408”+01 —00.0013 -—00.0030 
+24.0 +0.408 +00.0008 -—-00.9989 -—2.563'+02 —00.0009 —-00.0019 
+30.0 +0.120 +00.0003 +01.0003 —1.025”+03 =—-00.0004 —00.0008 
+36.0 —-0.178 —00.0003 —01.0004 =—4.101°+03 +00.0002 +00.0003 
+42.0 =0.461 —-00.0008 +00.9989 —1.64O“+OH +00.0007 +00.0014U 
+48,0 —-0,702 =—00.0013 =—01.0017 —6.562"+04U +00.0012 +00.C023 
+54.0 0.881 —-00.0916 +00.9278 —2.625”+05 +00.0916 +09.0031 
+60.0 —-0.980 -—-00.0018 -—-01.0024 =—1.050”+06 +00.0018 +00.0036 
+66.0 0.993 —00.0019 +00.9975 =4.200"+06 +00.0019 +00.0037 
+72.0 =0.917 —-00.0017 -—-01.0023 —1.680°+07 +00.0018 +00.0035 
+78.0 —-0.758 =—00.0015 +00.9981 -—-6.719“+07 +00.0015 +00.0031 
+84.0 —-0.533 —00.0010 —01.0014 —2.688°+08 +00.0011 +00.0023 
+90.0 —-0.259 —00.0005 +00.9993 —1.075”+09 +00.0006 +00.0013 
+96.0 +0.038 +00.0000 —-00.99939 =—4.300"+03 +00.0001 +00.0002 





SL, 


Generalizaties 


De voorafgaande beschouwingen werden gegeven voor een eerste 
orde differentiaalvergelijking, die nog lineair was ook. In 
principe gaat vrijwel alles net zo door wanneer f niet lineair 
is in x of wanneer in 61 x en f(t‚x) vektoren zijn inIR, d.w.z 
voor een stelsel van n eerste orde differentiaalvergelijkingen 
voor n onbekende furkties. Zoals bekend kan men een hogere orde 
differentiaalvergelijking altijd tot zo'n stelsel reduceren, 


en dat kan ook nog met stelsels hogere orde vergelijkingen. 


Bles 


Hoofdstuk 4 Wortels van vergelijkingen. 


Inleiding. 
In dit hoofdstuk beschouwen we het probleem een vergelijking 
f(x) = 0 op te lossen, waarin f een reële funktie van een reële 


veranderlijke is. Bekijken we bijv. 
Cl.) XÌ + 3X - U = 0 


Het linkerlid f van geze vergelijking 
is grafisch weergegeven in fig. 1.1. 
Een idee om deze vergelijking bena- 
derd op te lossen is het volgende: 

Zij x, in de buurt van de wortel. Bena- 


der dan de funktie f in het punt X, ZO 

goed mogelijk door een lineaire funktie 

g, (in de prent is dit gedaan voor Ks 

2); de grafiek van 8, Ìs dan de raaklijn 
aan de grafiek van f in het punt Cp sf xd). 
De vergelijking Ep = U is dan een he- 


naderde vergelijking, waarvan de wortel 





Xx, gemakkelijk te bepalen is 


Bo) = fx) + EE TTR Ae 


f(x) 
(1.2) X, ed Xo en Erk 


Men herhaalt nu het proces met X, Ì.p.v. X, enz. Zo krijgt men 


dus een getallenrij {x;} die voldoet aan 


Ex) 
CS) Xi+4 = BT Erle Le sr de Ar ten 
1 


Dit is het bekende proces van Newton(-Raphson) en uit de prent 


is het wel duidelijk dat dit proces heel snel convergeert, omdat 
in de buurt van de wortel de funktie g, (analoog aan g, gedefini- 


eerd voor x, i.p.v. Xp) een veel betere benadering voor f is 


1 
dan 8, en g, weer beter dan 8,» enz. 
Een proces als (1.3) waar men telkens m.b.v. een benadering 


Xx; een betere benadering Xx. bepaalt, heet een iteratief proces, 


1+1 


algemene gedaante 


CT) Kaa ò xs) 


Het is natuurlijk ook denkbaar dat men bijv. uit xj en Xi_4 


een nieuwe benadering x, bepaalt, of algemener 


1+1 


Cl Erg P DOK oger eere sj ) 


Processen van het type (1.4) noemt men 1-staps iteratieve processen; 
die van type (1.5) meerstaps iteratieve processen. Blijkbaar is 


Newton's methode een 1-staps methode met 


tn EED 
Lb HEX) = X raes) 
Men hoopt dat de rij {xj} naar een wortel a van de vergelijking 
f(x) = 0 convergeert enalsdefunktie & continu is in «a (in het 
geval (1.6) ziet men dat dit al zo is als f'(a)#0) dan geldt 
blijkbaar a = dla). M.a.w. elke eenstapsmethode transformeert een 
vergelijking f(x) = 0 naar de gedaante x = &(x) en lost deze 
vervolgens op met de methode van successieve substitutie 
Aang Te 
Deze opmerking leidt sommige schrijvers van elementaire boeken er 


toe een vergelijking als (1.1) maar op een ad-hoc manier naar de 


gedaante x = &(x) te transformeren en daarop successieve substitutie 


— 48 =— 


toe te passen. Bijvoorbeeld: 


zl 
(1.7) - X 2 E23 
en men krijgt dan beginnend met x, = 0, de rij 
(1.8) Ks 0, #, Ss Le3Is Ha 2 Urls Ry > dele Hj * Oe, 
Kg © 1.038, Xx, = O0 GB ee ae 


hetgeen een convergent karakter heeft, alhoewel langzaam. Maar 


wee de ongelukkige die eens geprobeerd zou hebben 





(1.9) xe F (h-x°) 

of 

10) xe 

Voor de laatste krijgt men, beginnend met x, = 1.01 (dus vlak 
bij de wortel) 

Edw Xg 5 1-01, xj 2 0.95, x, = 1.28, Xx, = 0.098, x, = 390. 


Men gebruike dergelijke ad-hoc methoden dus niet. 


52. Grafische interpretatie 
De methoden van successieve substitutie heeft de volgende 


grafische interpretatie 


53. 


3A 





vanuit 
naar y 
ant EL. 
en 2.4 


In het 


PLE. 28 fFis.dek 


X, Verticaal naar de grafiek van $, dan horizontaal 
= X»dan verticaal naar de grafiek enz. Men ziet dat er 
2.1 en 2.2 sprake is van convergentie, en in fig. 2.3 


van divergentie. 


volgende zullen we bekijken wanneer successieve substi- 


tutie convergeert, en zo ja, hoe snel dit gebeurt. 


Convergentie en convergentiesnelheid, 


Om inzicht te krijgen in het gedrag van iteratieve processen 


definiëren we 


{84 1) 





3B. 


het verschil tussen de i-deiterand en de wortel a, d.w.z. de 
fout in Xe Voor het proces van successieve substitutie 


geldt met de middelwaarde stelling van de differentiaalreke- 


ning 

iel dx) - dla) = P' CE) (x;-a) 

E tussen x; en a. hanen bla) = A; Xi44 7 dx) vinden we zo 
(3.3) Xaig T &P ò' CE) (x,-a) 7 
ofwel zn 
Gig * HOM, - 


Als we aannemen dat $' continu is in a, dan krijgen we voor h; 


klein, d.w.z. X; dichtbij a 


(3.5) hi,g  P'la)h, 


zodat bij elke stap de fout ongeveer met $'(a) wordt vermenig- 
vuldigd. We verwachten nu, dat voor |b'(a)| < 1 het proces con- 
vergeert als er dicht genoeg bij de wortel gestart wordt. Dit 
zullen we in stelling (3.7) bewijzen. Voor 0 # |é'(a)l <1 den 
noemen we $'(a) de (asymptotische) convergentiefactor van het 
proces. (Het geval b'(a) = 0 sluiten we uit omdat de mededeling: 
„bij elke stap wordt de fout ongeveer met 0 vermenigvuldigd! 
weinig zegt). 

Uit het voorgaande volgt dat het gunstig is als b'(a) klein is. 
Dat doet verwachten dat processen met b'(a) = 0 extreem snel 
convergeren. Zulke processen bestaan: ga na dat voor het Newton ais 


proces geldt &'(a) = 0 als f'(a) # 0. 


We zullen het voorgaande nu gaan formaliseren. 


Eerst een kleine generalisatie van (3.4) die een direkt gevolg 
van de formule van Taylor is. 


Stelling 3.6. Laat & k maal continu differentieerbaar zijn, met 


bla sun © sE ed = 0 en srad F 0. Dan is voor zekere 
funktie W met lim wt) = 0 
t->0 
_ gkd k 
hij 5 lb Ca) + plh‚)Ih;. 


Dit motiveert de volgende: 


Definitie 3.7. Een iteratief proces noemen we van de orde u 


als geldt 


[n 


= u 
gaal = IA + WD Inl 


AF 0, lim wt) = 0 
t>0 


Ga zelf na dat A en u door deze eisen eenduidig bepaald zijn. 

Een proces van de orde 1 resp. 2 noemt men ook wel een lineair 
resp. kwadratisch (convergent) proces, met voorbijgaan aan het 
feit dat er misschien helemaal geen convergentie optreedt. Het 
gedrag van lineaire en hogere orde processen verschilt aanzienlijk. 


We gaan dit nu onderzoeken. 


Voorbeeld Onderstaande tabel geeft wat resultaten voor het geval 


v(t) = 0. Maak zelf zo'n tabel voor A = 1, u= 3. 

















In stelling 3.9 zullen we zien, dat de tabel A = 1, u = 2 


representatief is voor alle kwadratisch convergente processen. 


De volgende stelling geeft aan wanneer een lineair proces con- 


vergeert: 


Stelling 3.8 Een lineair proces met |A| < 1 convergeert locaal. 


(d.w.z voor x, voldoende dicht bij a geldt lim x, = a). 
joo 


Bewijs: Voor ê klein genoeg is er een 0 < B < 1 zodat |A + w(t)l < B 


als 0 < |t S< 8. Laat dan nu I = (a-8,4+8). Voor x; € 1 geldt 


Eid - al < B Ix‚val <8 
dus ook Xi44 E I. Derhalve geldt : als Xx, É Il dan ook Xi; E I 
voor elke i en |x; - al < B [x‚- af. Omdat 0 <B <1 is 
lim |x, - al = 0 d.w.z. lim x; = a. 
Loo ie hed: 5 


Laat nu zelf zien dat voor |A] > 1 een lineair proces dat dichtbij 
de wortel gestart wordt, de neiging 
heeft om van de wortel weg te lopen. 
Eenmaal ver van de wortel afgeraakt, 
kunnen andere effecten een rol gaan 
spelen, waardoor we toch weer bij 

de wortel terug kunnen komen: 

in nevenstaande figuur is toevallig 
X, = A, hoewel in het begin 


[xd „al > [x, = al omdat 





|A] = |é'(a)| > 1. In de numerieke praktijk zal door afrond- 
fouten x; toch wel weer uit a weglopen, zodat dit voorbeeld 


eigenlijk niet eerlijk is. 





Shen 


Iteratieve processen met u > 1, zgn. superlineaire processen, 


gedragen zich netter 


Stelling 3.9. Een superlineair proces convergeert altijd; ongeacht 
de waarde van A. 
Bewijs : Laat |A + W(t)| S B voor |t| < 8.Door de schaal op de, 


x-as te wijzigen blijkt er een prettige relatie tussen de fout 


in stap ì en die in stap i+1 te voorschijn te komen : beschouw nl. 
A. 

koe BPS, 

1 zl 

_i_ 
Met 8' = BĲ Î & geldt dus voor k, € I' = (-8',8') dat 

< u 
grat SS Hel 


1 zÂ N Li 
Voor k, € I' en [k, | < 1 is dus ook k;,, € I' en [ks 44! Ede ALE 


i 
k, € I' en |k,l < 1 geldt dus [k; | < ie, 15 . Omdat u° >» voor i > @ 


geldt dat lim k, = 0 en dus lim h‚ = 0. Voor h, voldoende klein 
Loo 1 jo L 
treedt dus convergentie op. 


Opmerking 1. We zagen reeds dat k; de fout is in de i-de iterand 
in een andere schaal op de x-as. De relatie Ik; 44! S lik, [* be- 
tekent dan dat voor [k, | < 1 het aantal goede decimalen achter de 
komma per slag met een factor u toeneemt. Zo zien we dat door toe- 
voeging van de schaalfactor uit het bewijs van (3.8) de zojuist 
gegeven tabel representatief wordt voor alle tweede orde processen. 
2. Onze verwachting dat een proces met b'(a) = 0 snel 
zal convergeren wordt door (3.7) bewaarheid : superlineaire pro- 


cessen convergeren dichtbij de wortel zeer snel. 


We willen twee superlineaire processen met elkaar vergelijken om 


te zien welk sneller convergeert. Wat we daar precies mee be- 


doelen zegt de volgende 


Definitie 3.10. Een proces P (met “fouten! h‚) heet sneller convergent 
dan een proces P' (met "fouten" hi) als voor elke h, en hj waar- 
voor de beide processen convergeren geldt dat h; S hi voor allè 


voldoend grote i. 


Een sneller convergent proces haalt een langzamer convergent 
proces dus altijd in, waar ook gestart, mits ze maar convergeren. 
Laat nu een proces P van orde u > 1 gegeven zijn, en een proces 
P'van orde u' > 1. Laten we aannemen dat u > u' en dat beide pro- 
cessen convergeren. Als we met proces P voldoende ver gevorderd 
zijn, dan winnen we -op de juiste schaal bekeken- per slag u maal 
zoveel goede cijfers achter de komma. VoorP!' is dit =op een wellicht 
iets andere schaal- u' maal zo veel. Men verwacht nu, dat het 
proces P sneller convergent zal zijn dan op het proces P!' omdat 
bij P het aantal goede cijfers sneller groeit dan bij P'. Enig 
nadenken leert dat dit niet verstoord wordt door het verschil in 
schaal. Voor de liefhebbers geven we een formeel bewijs van een 


en ander. De stelling luidt 


Stelling 3.11. Laten twee processen P,P' gegeven zijn van de orde 

u‚u'. Neem aan dat u > u'. Dan is P sneller convergent dan P!, 
OR: $ t 2 t ks 

Bewijs: Laat P fouten h; hebben, en P' fouten hi . Wegens h, > 0 


is er volgens het bewijs van stelling 3.9 een rij k; = ch, zodat 


vanaf zeker rangnummer i, 


< u 
geldt Ik, | < 1 en [ks 44! [| 5 
Analoog aan het bewijs van 3.9 toont men echter aan dat er een 


symmetrische omgeving J van OQ en een getal c!' bestaan zodat voor 


SU, 





_- 55 =— 


t 
de pij k! = c'h! geldt |k!‚,| > kt {P mits k! € J. Omdat ook 
aL sb 1+1 1 1 


h! > 0,liggen alle k! vanaf zeker rangnummer Ln in J, en geldt 
1 d 





t 
dus vanaf dat rangnummer Des, 4 | AE . Zij N = max (i,sig) 
Dan geldt voor i > N en voor B zodanig dat [kil = [iel 
„(d-N) 
De; l  Ikyl EN wen gt Ci-N) 
SS q 

t CIN) N 
De | H 

B 
en hierin geldt pi gur GD > wo voor i *@v wegens u > u'. 


Dus k;/k; > 0 wegens |k‚| < 1, dus h;/hì * 0 wegens h;/h; = (e'k)/(ek;). 


N 


Newton Raphson 
Zoals we zagen is het gewenst een vergelijking f(x) = 0 te 
transformeren naar de gedaante x = b(x);, en de meest voor de 


hand liggende wijze is dat men neemt g(x) = Xx t f(x): Maar dan 


convergeert de rij Xx. 
8 J Xi+1 


waarvoor geldt f'(a) > 0 of f'(a) S - 2. 


= dx) niet meer naar een wortel a 


Vandaar dat we proberen dx) = Xx + Alx)f(x), Ax) geschikt te 
kiezen. Dan is onder aanname van voldoende differentieerbaarheid 
b'(a) = 1 + XAla)lf'(la), zodat A voldoet als Ala) = -1/f'(a). Hier- 
aan wordt bijv. voldaan als A(x) = =1/f'(x) en men krijgt 
Newton-Raphson terug (zie (1.3)). 

We zagen al dat de orde van NR groter dan 1 is als f'(a) # 0. 

Om de orde nader te bepalen moet men hogere afgeleiden van & uit- 
rekenen in a. Zo is &"'(a) = £f"(a)/f'(a), dus voor f''(a) £ O0 

is de orde precies 2. Bij het uitrekenen van &' verschijnt ook 
een term £f'''., Dit introduceert echter een extra differentieer- 


baarheidseis voor f. Hoewel men daar in de numerieke wiskunde niet 


zo zwaar aan tilt, geven we toch ook een berekening van de orde 
van NR waarbij geen hogere afgeleiden dan de tweede voorkomen. 
Tevens geeft deze afleiding een precieser verband tussen A a 
en Xx; — a. De afleiding berust op een merkwaardig gebruik van 
de Taylor reeks, nl. niet met als steunpunt a (wat men zou ver- 
wachten) maar Xi; Onder de aanname dat f'(x) bestaat op een 
omgeving van a geldt 


(a-x.)? 


(u.1) 0 = fla) = he er ant di 5 


Voorts (1.3) 


(u. 2) 0 EOD) grid fl Ox) 
Aftrekken geeft 

8 £"CE) 
(u. 3) Xj44T 0 = (xjva) EA EN 


f(a) _ 
2E! la)? HK 7 Ô 


en dit is inderdaad van de gedaante (3.7) met A 


als f'(a) # 0 en f(a) # 0 en f" continu in a. 


We zien hieruit dat f(a) = 0 een gunstige invloed op de con- 
vergentiesnelheid zal hebben terwijl f'(a) = 0 het omgekeerde 
effekt zal hebben. Inderdaad gaat het gunstige convergentie 
karakter van NR verloren als f'(a) = 0. Dit is onmiddelijk 


duidelijk uit het volgende voorbeeld 


k 
Xs 
___k k a Aen ziet 
f(x) = Xx. Dan is Xia4 5 *i me EA (1 KX zodat men ziet 
1 


dat bij deze k-voudige wortel er per stap slechts een fraktie 


rl Ke 


(3) 


van de fout af gaat. Men kan aantonen dat NR algemeen ook locaal 
convergeert in het geval van een meervoudige wortel. De con- 


vergentie is dan lineair. Bewijs dat zelf. Door a in 


- n f(x) 
dx) = x=-a FTC) 





1 1 
mf (EN), Ela ,x;) of Ge, 


s OL 


k 


85. 


SA 


passend te kiezen, krijgt men dan wel een kwadratisch conver- 


gent proces. Ga dit zelf na, als f een k-voudig nulpunt heeft. 


Tenslotte keren we terug naar het voorbeeld in 81, Newton- 
Raphson toegepast op f(x) = x°+3x-4, startend in Xg = 2, waarvan 


we de resultaten geven in de volgende tabel: 


DT irr 























: 2 
4 Xi; h;/(h;_,) 

jo 2.0000000000000 = 
1 1.3333333333333 0,33 
2 1.0488888888889 O.bb 
3 1.0011751554604 0.ug « 
hj 1.00000069C2246 0.51 
9 1.0000000000002 O.u2 

- 1.0000000000000 = 

en we zien de fraaie overeenstemming van h;/hi_4 wet Ela) = Db, 


Koorden Newton 


Voor eenvoudige funktiesf is NR een zeer handzaam proces. Als f 
echter een ingewikkelde funktie is pleegt f!' nog veel ingewikkel- 
der te zijn, met overeenkomstige stijging in rekentijd. De vraag 
doet zich dan voor of het nu echt wel nodig is de funktie f£ in 


het punt x, zo goed mogelijk door een lineaire funktie g, te be- 


0 
naderen en of het niet wat minder kan; tenslotte geeft g, ook niet 
de exakte wortel. Het lineair Ë 
interpolatie polynoom in Xx, en 


een naburig punt x_, lijkt ook 


n 
niet zo'n gekke keus. Men kiest 
weer x, als nulpunt van g, en 


5, als lineair änterpolatdie 


polynoom in Xx, en X,» enz. 





Zie fig. 3.1. In feite betekent dit dat men in het NR proces 


£'(x,) vervangt door 





Slad en Phed 
DL Med 
Men krijgt dus 
BE Od 
As ek IE ETE Es Cet ahd 
1 1-1 


Dit proces noemt men koorden-Newton (ook wel secant methode of 


regula falsi). Het proces heeft de vorm Xi44 = dx; s Xi) en 1s 


dus een 2-staps methode. Hiervoor kan men het convergentiegedrag 


niet meer zo bepalen als bij NR. 


Onder de aannamen f'(a) # O0 en f' bestaat op een omgeving van 


‚a geldt (zie (4t.1) en (4.2) van hoofdstuk 1) 


(5.2) O0 = fla) = f(x) + Ca-x; df (xy sj) + Ca-x;dlarx; fx; xj gp) 


Verder (zie (5.1)) 


(5.3) Oz fx) + Os = Xi df (xs oxsg) 


1 1 


aftrekken van beide formules geeft 


| aon ns 8 ETE) 
(5.8) Xi4g 7 AF Ca-x;) (a Xie) ZET Cn) 
Ofwel 
EEE) 
ke hisa * Tin) * Piter: 


(zie de gelijkenis met (4H.3)). 


Berekening van de orde van KN is niet eenvoudig. Voor liefhebbers 


wordt in het volgende een afleiding gegeven. Wij vermelden hier 


as 


slechts dat de orde van KN gelijk is aan A = Za + /5) = 1.62. 


In feite is 


h. 


A-1 
RL 


E en ‚ 
(5.6) In; 44! = |A | 


SB, 





Opmerking. Uit 3.5 ziet men dat voor fatsoenlijke funkties een 
Éênstapsproces altijd een gehele orde zal hebben. Voor meerstaps- 
processen hoeft dit blijkbaar niet te gelden, zoals we zojuist 


zagen. 


Orde KN voor de liefhebbers. 


" 
1) Convergentie Zij Kaar < B op een symmetrische omgeving 


En LB Ä 4 < h!'h! 
van a. Dan geldt voor h; | Bh, | dat h; | Bh, | dat hi,4.S hihi, 
als x; en Kid in die omgeving liggen. Als x_q en X, nu ook nog zó 
dicht bij a liggen dat h, en ho <1 dan is de rij hjshishzer: 
monotoon dalend, heeft dus een limiet, en die kan niet # 0 zijn. 


Dus er is locale convergentie 


2) Heuristiek over het gedrag van {hj}. Zij f'" continu en # 0 


p or AR n 
in a. Zij A = DEL) Dan geldt dus h; 44 = (1+b JA hijh;_4 met 


ò; > 0 wegens de zojuist vastgestelde convergentie. Dus, met 


k. 
re 


|Ah; |: 
si 


(5.7) kos 7 Mtb) kijks;4 


We vermoeden dat de rij (u; } gedefinieerd door 


(5.8) Ussg * UjUsoge Uig k_,» us = k, 


een indruk van het gedrag van {k,} geeft. Logaritme nemen geeft 


Er SMA = 0 met LE log u; Vv s JOB Ke 


rt a Lel e logk_4» hd 


a o 


Deze betrekking voor {vj} is een zgn. drieterms recurrente be- 
trekking en alle rijen {v;} die eraan voldoen laten zieh eenvoudig 


Az 


aangeven met behulp van de wortels À = 5/5) = 1.6 en 


uz Atte) z =0.6 van de zgn. karakteristieke vergelijking 


2 
x?-x-1 = 0. Voor deze wortels geldt nl. À 
à 


mi me 8 ; 
en u Ti u = 0, zodat elke Pij {v‚} met Vi = art+gur, a en B 


it1_jdjiel : KE bene) 


Ht 
Lan) 


willekeurige konstanten, aan de recurrente betrekking voldoet. 
Omgekeerd is echter elke rij {v;} die aan de recurrente betrekking 
voldoet op deze wijze te schrijven, want neem maar a en & zo dat 
atB = vo » aA+Bu = v, en overweeg dat een rij die aan de betrekking 
voldoet volledig bepaald is door twee opeenvolgende elementen eet 


ervan. 
E 1 sl ' n _ i-1 . 
Er geldt dus Vi © ah” +Bu. Dan is virÂvi.g z Bu (u-A) en dit 
U. 4 


gaat naar OQO voor Ì > @ wegens u=-0.6. Dus lim log EE kr 0, 
10% u. il 
Aes 


À 
4 We vermoeden dus ook dat k; = (1+o(1))k;_ 


en dat zou betekenen dat KN de orde A > 1.6 heeft 





zodat u. = (1+to(1)) Mis 
â sf 


d 
K; 
3) Precies gedrag der k; (en dus h‚). Beschouw w; = (waarvan 
k. 
gt 
we dus hopen w; > 1). Dan volgt uit (5.7) (bedenk A(1-A) 


„1 
en u = 1-À) 


1d 1-A, -ACA-A) 1-À u 
5 Kgs Cdi “kies CAtpjduy D= (Arde n 


W 
) 


das (ltb‚)k 


weer log nemen geeft met in log wi _ 


Piat log C+) + HY à v; + Hrs 


i hs 
met lim p, = 0, Dus (met volledige inductie) 
s 1 
(5.9) Baar © vs; + UP; _4 F zen PMD, 


Zij nu e > 0 gegeven en zij lv; 1</3 voor i > i,. We splitsen 
de som in (5.9) bij de index i, en schatten de beide delen we 


afzonderlijk af. Dan geldt 


oC 


en 61 ee 
Let, ietgt dt Ld 
rz) S Zaldelnlde ttl l + lul [ielul+...+lul Lmax hd | 
J 
5 [ed Lul max bw; [1/A-lulD GE 
voor Ì groot genoeg wegens 1-|u| = 0.4. Dus lim Ek 0 dus 
lim wss 1. Hiermee zijn we klaar. 


Opmerking. Bovenstaande methode om alle oplossingen van een re- 
currente betrekking voor te stellen m.b.v. machten van de wortels 
van een karakteristieke veelterm is algemener toepasbaar. 


Voor details zie Numerieke Analyse I en III. 


Getal lenvoorbeeld 
We besluiten deze paragraaf met KN toegepast op het voorbeeld 


uit sle Uit (5,6) volgt det 


lim De 
100 


10 en 
log h;/ log hijs * A 


In de tabel zien we dat dit quctiënt aardig tot A nadert. 










10 10 
log h;/ log hi 








„0000000000000 
.obSHSHSUSU5US 
„1986754966887 
„0467569961944 
„OOUUU491206U1 
„0001029408471 
„0000002284535 
„0000000000118 
„0000000000000 
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Koorden Newton toegepast op f(x) = x°+3x-4, startend met Ha 8 





86. 


BT. 


Vergelijking van NR en KN. 


Eigenlijk bevat de vorige paragraaf een teleurstelling omdat 
de orde van KN lager is dan die van NR, wat men eigenlijk niet 
verwacht zou hebben, aangezien voor Ì > @ de punten x; steeds 
dichter bij elkaar komen te liggen zodat de koorde op den duur 


zeer goed de richting van de raaklijn aanneemt. 


NR convergeert dus sneller dan KN (stelling 3.11). De lezer zij 
echter op zijn hoede : dit wil zeggen dat op den duur NR met 
hetzelfde aantal stappen dichter bij de wortel komt dan KN. 

Het betekent echter niet dat NR ook met een zelfde hoeveelheid 
rekentijd op den duur dichter bij de wortel komt dan KN, want 
KN heeft minder tijd per stap nodig. Omdat f' meestal veel in- 
gewikkelder is dan f is het heel redelijk te zeggen dat KN per 
stap doorgaans hoogstens de helft van de rekentijd nodig zal 
hebben die NR per stap nodig heeft, oftewel dat KN in de tijd 


van één stap NR twee stappen doet. 


Nu geldt voor KN op den duur h. a Chí, \ = 1.6 (zie 85), dus 


11 
hi As gard „A ‚ hetgeen betekent dat wanneer we een proces P 
definiëren waarvan één stap bestaat uit twee stappen KN, dit 
proces P van de orde A? = 2.6 is. Dit proces komt dus op den 
duur met evenveel stappen als NR dichter bij de wortel, en heeft 


per stap hoogstens evenveel tijd nodig als NR. Dus in rekentijd 


gerekend is KN sneller convergent dan NR. 


Globale convergentie 


Tot dusverre hebben we slechts het gedrag van methoden vlak bij 


de wortel bekeken (het zgn. asymptotisch gedrag). We willen nu 





iets zeggen over water vervan de wortel gebeurt. Allereerst 


merken we op dat uit plaatjes duidelijk is dat NR en KN monotoon 


Xo X0 


convergeren als de startpunten liggen in een interval met de wortel 
als eindpunt waarop de zraftek van f de bolle zijde naar de 

x — as keert. Dit is ook exact te bewijzen; de aanname luidt dan 
dat f'"'(x)f(x) #» 0 op zo'n interval. 

Dit zegt echter nog niets over de snelheid van het proces 

ver van de wortel. 

Voorbeeld 1. We bekijken de vergelijking x° - a = 0. De NR itera- 


tieformule is dan te schrijven als 


ie a 1 
Ted Raa 5 Cx + ri 12 (ga na!) 
Neem a = 1 en x, = 10° (ga na dat monotone convergentie naar de 


0 
wortel + 1 gegarandeerd is). Zolang nu de iterand groter is dan 5 
(d.w.z. Ìi < 15) is er een vrijwel lineair verband tussen h, en 


h. 
h Ae! ne ee, zodat NR in het begin niet sneller convergeert 


zal … 1+1 2 
dan een halveringsproces. (Een halveringsproces ter oplossing 
van f(x) = 0 is een proces van het volgende type : laat f continu 
zijn en tussen a en b één wortel a hebben, terwijl f(a) en f(b) 
verschillend teken hebben. We kiezen dan ec = (atb). Door naar 


het teken van f(c) te kijken, zien we of a in (a,c) of in (ec;b) 


ligt (of azc). De wortel a is nu in een 2 maal zo klein interval 


88. 


8A 


gelocaliseerd). 
Als er echter dicht bij de wortel wordt gestart zien we duidelijk 
het kwadratische convergentiekarakter. Met x, = 1.1 vinden we 


achtereenvolgens hj = 101 ; h, = 4.5 Xx 10%; h, 2 10° en 


11 


h, = 5 Xx 10 (vergelijk dit met het theoretische verband tussen 


hi44g en h‚) en we zien, dat nu iedere NR slag een verdubbeling. 
van het aantal goede decimalen geeft. 


2. Voor f(x) = (x?-1)(x?-4) en Xo 5 10° convergeert 


NR aanvankelijk lineair met factor à. Ga dit zelf na. 


3. Laat f(x) = e“-1. Door X, groot te nemen ziet men 
dat NR aanvankelijk zelfs willekeurig langzaam kan convergeren.- 
Hetzelfde geldt wanneer f een asymptoot heeft en men start vol- 


doende dicht bij deze asymptoot. 


Men ziet dus dat men op grote afstand van de wortel weinig plezier 
heeft van de hoge orde van een proces. Dat plezier begint pas wanneer, 
in de juiste schaal gemeten, h; duidelijk kleiner dan 1 is ge- 
worden (stelling 3.9). 

Vandaar dan ook dat men vaak met ruwe, goedkope methoden een be- 
trekkelijk klein interval bepaalt waar de wortel in zit, al- 

vorens deze met een lokaal snel konvergente methode echt goed te 
gaan benaderen. Zo'n ruwe methode is bv. het eerder genoemde 


halveringsproces. 


Gevolgen van afrondfouten 


Iemand die met de hand rekent, rondt gewoonlijk op een ad hoc 


manier af. Wij zullen aannemen dat dit systematisch gebeurt, en wel 


als volgt: kies een natuurlijk getal M‚, en rond elk getal af op 
het dichtsbijzijnde getal m. 10P, m en p geheel, [ml < 101 
(getallen m.10P met [ml < 10% hoeven dus niet meer afgerond te 
worden). Dit geldt zowel voor de begingetallen als voor het 
resultaat van elke optelling, aftrekking,vermenigvuldiging en 
deling. We noemen dit : rekenen in M decimalen nauwkeurigheid. 
Deze wijze van rekenen komt overeen met die van een computer 

(zij het dan dat die pleegt af te ronden op getallen m.2P). 

Door op deze wijze het gebeuren in een computer te simuleren, 
hopen we inzicht te krijgen in de gevolgen van het rekenen met 
eindig veel cijfers bij het bepalen van wortels van vergelij= 
kingen. 

Het rekenen in eindig veel cijfers betekent dat bij iedere hande- 
ling waarbij een getal optreedt dat niet van de gewenste gedaante 
is (zelfs al bij het opslaan van zo'n getal in de computer) 
afrondfouten gemaakt worden. Men kan eenvoudig inzien dat de 
relatieve afrondfout in de grootte orde van oink is, wanneer M = 12. 


Dit laatste zullen we in het volgende steeds aannemen. 


Bij het numeriek rekenen is men gedoemd afrondfouten te maken. 
Dit betekent dat men de funktiewaarden slechts met een zekere 
onzekerheid te pakken krijgt. Stel dat deze onzekerheid e is, 
dan impliceert dit een 
onzekerheid in de nulwaarden. 
In fig. 8.1 ziet men dat nu 
alle getallen än het Mmrter= 
val [B‚ylmogelijke nulwaar- 


den zijn geworden. We noemen 





daarom [@&‚yl het onzekerheids=- 


interval voor de nulwaarde. 





8B 


Een andere wijze om dit zelfde 
interval te karakteriseren ziet 
men in fig. 8.2 : het is het in- 
terval waarop geldt |f(x)| Se. 
Uit fig. 8.2 zien we dat in 
eerste benadering geldt 


Yy-aza=-gTe/rlf'(oa)|. mits 





£'(a) # 0. Als f'(a) = 0 maar 


1 
2 


f'(a) # 0 dan geldt met Taylor f(x) = F"(E) (x-a)?, zodat nu 


IEGOL ze y-aTa=Bg = (IETTEMAT (zie fig. 5.3), waarbij 
men moet bedenken dat /e voor kleine e veel groter is dan € 
zelf. 

We bekijken eens wat dit laatste betekent voor de eenvoudige 
funktie x? - 2ax + a?, die voor x = a een dubbele wortel heeft. 
Voor x in de buurt van a zijn x°, ax en a? ongeveer even groot. 
Dit betekent dat e = ua? 10 1 (vgl. 8A). Voorts is f''(a) = 2, 
zodat y - az gado waaruit men ziet dat alhoewel men slechts 
relatieve afrondfouten maakte van 10, de onzekerheid bij de 


wortel 6a. 10° geworden is, d.w.z. dat men 6 cijfers verloren 


heeft. 


Deze onzekerheidsintervallen zijn ook van betekenis als men 

de nulwaarden numeriek wil benaderen. Immers, voor vrijwel alle 
funkties (bijv. kwadratische funkties vormen een uitzondering) 
zijn alleen maar iteratieve methoden bekend om nulwaarden te be- 
naderen en deze leiden in principe alle uit het teken van de 
berekende funktiewaarde af of de nulwaarde links of rechts ligt 
van het punt waar men zit (zie bijv. NR: als f'(x) > 0 en 

fx) > 0 dan is X; S Xi): Dit betekent dan echter dat men de 


+1 


verkeerde kant uitgestuurd wordt als de berekende funktiewaarde 


8C 


8D 


het verkeerde teken heeft. Derhalve mag men van geen van deze 
methoden verwachten dat ze iets beters doen dan willekeurig dicht 
in de buurt van het interval [B‚yl te komen. Ook ziet men dat 
wanneer de methode toch in het interval komt, men moet verwachten 
dat een vreemd heen en weer zwalken van de opgeleverde puntenrij 
optreedt. Wanneer we weten dat de rij punten van NR monotoon. 
naar de wortel gaat kunnen we dit verschijnsel op de volgende 
wijze uitbuiten: zo gauw de op de computer verkregen rij op- 
houdt met stijgen (of dalen) stoppen we de berekening. De laatst 
verkregen waarde ligt (hopen we) redelijk in de buurt van het 
onzekerheidsinterval. Dit blijkt in de praktijk goed ús werken. 


Op deze wijze is het NR proces in het voorbeeld in $4 gestopt. 


Op nog geheel andere wijze maakt men met afrondfouten kennis bij 
KN. Immers, als men vlak bij de wortel is aangekomen kan de waarde 
van fx) je fg) geheel door afrondfouten worden bepaald, en 
daarmee wordt de richtingscoëfficiënt van de koorde volkomen 
onbepaald. Deze zou bijv. zelfs 
wel O0 kunnen zijn, waarmee de 
zaak geheel uit de hand loopt 
(zie fig. 8.4). Men moet bij 

KN dus wel van ophouden weten. 


Gezien de snelle toename van 





het aantal goede cijfers vlak 
bij de wortel (zie 585) is een vuistregel dat men ophoudt zodra 


IE < 10 P [x; | waarin p zoiets is als 5 van het aantal 


ij Xi 3 


cijfers waarmee gerekend wordt (zie BA : p = 5m). 


Als laatste zeggen we nog iets over wortels van veeltermen, aan= 


gezien men deze in de praktijk vaak ontmoet. Het ligt dan voor de 


hand om, als men eenmaal een wortel a heeft gevonden, deze uit te 
delen en verder te gaan met de aldus verkregen (n-1)-ste graad 
veelterm g (g(x) = f(x)/(x-a)) (minder werk en meervoudige wortels 
komen er vanzelf ook meervoudig uit). Echter, als men een 
polynoom deelt door Xx - a, waarin a een benaderde waarde is voor 
een wortel a van dat polynoom, zal het gereduceerde polynoom 
t.g.v. afrondfouten bij het uitdelen en de fout in a afwijken 

van het exakte polynoom f(x)/(x-a). Belangrijk is nu, in hoeverre 
de wortels van het gereduceerde polynoom hierdoor verknoeid zijn. 
Men kan aantonen, dat het op de gebruikelijke wijze uitdelen 

van de wortel met grootste modulus ernstige fouten kan geven 

en dat dit bij uitdelen van de wortel met kleinste modulus niet 


het geval is. We laten dit zien aan een 


Voorbeeld Het polynoom x°-999x?-1001x-1 heeft als wortels 


z= =-1, a, = -0.000999.... Rekenend met 


a, = 1000 .001....., a p 


1 2 


getallen van de gedaante m.10P,|m| < 10°, m en p geheel, is a = 1000 


de allerbeste waarde die er voor a, bestaat. Uitdelen 


x-1000/x% -— 999x? - 1001x == I\x? + x-1 
x3 -1000x? 
x2 
x2 = 1000x 


=X 


en natuurlijk gaat het niet op maar dat mocht men ook niet ver- 


wachten. Evenwel heeft de na deling verkregen veelterm x° + x= 1 
de wortels = 5 + 5 /5 = 0.6 en - 1.6 en dat is niet zo goed. Had 


men evenwel alsvolgt uitgedeeld 


„_ 69 =— 


1000-x/-1 - 1001x - 999x2 + XIN -0001=1.D0ix-x? 
-1 +0.001x 


… 400 Te 


= 1001x+1.001x? 
te ien lp 
= 1000x? 


dan heeft de na deling verkregen veelterm x? + 1.001x + 0.001 de 
wortels = 1 en - 0.0001 en dat is heel goed. Maar deze wijze van 


uitdelen komt er dan ook juist op neer dat men = uitdeelt uit 


=xS — 1001x? - IIR F 1, dus het polynoom met wortels de, a en Ee 
1 2 


en hierin is En juist de in modulus kleinste. 
1 
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Stel dat het getal a een benadering is voor een op een of andere 
wijze gedefinieerde grootheid a, dan noemen we de fout in a: a - a. 
De absolute fout in a zal zijn : |a-al. Onder de relatieve fout in 
a verstaan we : En : 

Vaak geven we de relatieve fout ook in procenten aan door uit 


te rekenen : relatieve fout Xx 100%. 


Foutfiliosofie : het percentage van de relatieve fout hoef je niet 


in grote procentuele nauwkeurigheid te kennen. Je bent er wel in 
geïnteresseerd of een fout 1% of 3% is, de mededeling dat een 
fout 1.1% is echter nauwelijks inhoudsrijker dan dat hij 1% is. 
Zo zullen we rustig zeggen dat het getal 0.2347, ontleend aan 
sinustafel in 4 decimalen, een relatieve fout van omstreeks 0.023 


heeft, alhoewel dit in feite 0.02132...% is. 


(a) Ga na dat de uitspraak "a heeft een relatieve fout e" gelijk- 


waardig is met a = a(1+e!), waarin |e'| = e. 


(b) Ga na als relatieve fout (a) = e‚ en relatieve fout (b) = n 
met €, en €, <1 
(1) relatieve fout (ab) S ete, 
K24 pelatieve fout (1/b) * Ee, 
(3) relatieve fout (a/b) S Ee +E, 
(ht) relatieve fout (atb) S max(e,;e,) als a en b hetzelfde 
teken hebben 


Bite 
(5) relatieve fout (atb) & 3 als a en b voldoen aan 





… 


2 z= -1+ô, ê > 0 5; dus bij het optellen van getallen die 
bijna elkaars tegengestelde zijn (d.w.z. 8 < 1) kunnen 


de relatieve fouten erg opgeblazen worden. 


(e) Ga na dat in (1),(3),(4) en (5) de relatieve fout best kleiner 


kan zijn dan de in het rechterlid genoemde grootheid. 





dk 


(d)* Ga na dat niettemin (1) en (3) realistische foutschattingen 
zijn in de volgende zin : voor elk tweetal getallen a en 8 
en elk tweetal niet negatieve getallen e‚ en e‚ bestaan er 
benaderingen a en b van a en B met relatieve fouten e‚ en €, 
zodat de relatieve fout in ab resp. a/b ongeveer gelijk is 
aan de waarde in het rechterlid van (1) resp. (3). 


Ga na dat dit bij (5) ook zo is mits ê klein. 


Hoe zit dit bij (4)? 


In deze opgave willen we nagaan hoe het staat met de nauwkeurig- 
heid waarmee men een functie kan benaderen als men over 4 waarden 


van een functie en zijn afgeleiden beschikt. 


Zij f Ux continu differentieerbaar. We vergelijken nu de uit- 

komsten op het interval [O,h] van 

(1) k-punts Lagrange-interpolatieop de steunpunten =-h, O‚h en 2h. 

(2) 2-punts Hermite-interpolatie op de steunpunten Os;h (d.w.z. 
gebruik f(0), f'(0), f(h) en £f'(h)). 

(3) Een Taylorreeks rond x = 0, afgekapt na de 3e graadsterm. 


Jh, afgekapt na de 3e graadsterm. 


(tk) Een Taylorreeks rond x 
a. Ga na dat in elk der 4 genoemde methoden de fout gegeven 
wordt door R; (x) an: m0) 
met TG) een be graads polynoom behorend bij methode (i). 
b, Laat zien waar 1; op [O‚h] extremen aanneemt (inwendig- of 
rand-extremen). Bereken de extremen. 


ec. Laat zien dat voor x € (O‚h) : mx) > m,(x) > 0. 


Schets nu a,s7,stTg en Ty op het interval [O‚hl. 


3). 


Ca 


Stel iemand wil een tabel van goniometrische functies opsteller, 


-4 


waarbij de functiewaarden hoogstens een relatieve fout van 3.10 


Rand 


bezitten, om daarin prettig lineair te kunnen interpoleren. 


a). 


b). 


ea). 


Met welke stapgrootte h moet men sin(x) op[0,7/u} tabuleren op- 
dat de absolute fout bij lineair interpoleren kleiner is dan 
10 *e 

De fout bij lineaire interpolatie van de sinus op de punten 


XX + h wordt gegeven door: 


_ sin(&) 


5 (x= Xo) Cx EK Slide 5 Shan + DJ 


0 


We benaderen nu sin(&) in deze formule door de zojuist met 
lineaire interpolatie verkregen waarde q van sin(x), zodat 


de fout nu benaderd wordt door 
== A/È (5 Kdl e= xp > hJ 


Bij het resultaat q, verkregen met lineaire interpolatie, tellen 
we nog de zo verkregen schatting van de fout op. 

We keijgen dus 5 AIT — blm == HJR =S HR =S HJJs 

Laat nu zien dat het eindresultaat behept is meteen fout van 

ten hoogste h3/g voor x € [0,7/yl. 

Ga tenslotte na hoe groot h moet zijn opdat lineaire inter- 
polatie + bovenstaande trucage een fout kleiner dan sn op= 


levert? 


Benader nu met lineaire interpolatie en de methode uit b) 
sin (0.45) gebruikmakend van sin (0.4) = 0.3894 183423 
sin (0.5) = 0O.4794 255386 


N.B. : sin (O.45)= 0O.4349 655341 


U). 


Lineaire interpolatie gaat uiteraard alleen goed als de functie 
zich goed door een lineaire functie laat benaderen. 

Laat nu tan(1.566) = 208.5 en tan(1.568) = 357.6 beide afgerond 
in hk cijfers gegeven zijn. Doel is nu tan(1.567)(= 263.4) in 
een vergelijkbare nauwkeurigheid te bepalen. 


N.B. : 7/3 = 1.5708). 


a). Men berekent tan(1.567) met behulp van lineaire interpolatie 
op de punten 1.566 en 1.568. Geef een redelijke ondergrens 
voor de fout met behulp van de theorie, 

Hoeveel goede cijfers kunt U dus verwachten in de zo bereken- 


de waarde van tan(1.567) ? 


b). We konden ook niet verwachten dat we goede resultaten zouden 
krijgen, want tan(x) = HET: en dit is bij de asymptoot niet 
goed met een lineaire functie te benaderen. Ons idee zou 
dus kunnen zijn : beschouw 1/tan(x) = cot(x) want die ge- 
draagt zich ongeveer als 7/9 - x en dat is wèl lineair. We 
verwachten dus via interpolatie op de cot(x}) en dan terug naar 


tan(x) betere resultaten. 


e). Benader nu cot(1.567) met lineaire interpolatie als cot(1.566) 
en cot(1.568) bekend zijn. Hoe groot is de relatieve fout 
(geef een bovengrens)? Dit is dan tevens de relatieve fout 


in de hieruit verkregen waarde van tan(1.567) 


(zie opgave 1) 





cot(x) Eostx) 
0.004796 





tan{(x) 









208,5 








0.003796 





263,4 









0.002796 





357,6 





Be 


B), 


We beschouwen een equidistante tabel met Ax = h voor een 


fe CGC Ta,bl, aen b edndig: 

a). Laat zien dat bij k punts interpolatie (k = 2 of 3) een abs. 
fout S el”, gemaakt wordt, en dat als h voldoende klein is 
deze grens realistisch is, d.w.z. dat de abs. fout de waarde 
en* (bijna) kan bereiken. 

b). Laat zien dat voor h > 0 de foutschatting bij 3 punts inter- 


polatie sneller naar O gaat dan die bij 2 punts interpolatie. 


Een tabel CTEELSTRE, van een functie f zou men ook kunnen 


Zed eat 


gebruiken als tabel ACTRESS REDEN ‚ hetgeen geïnterpreteerd kan 


caf 
worden als een tabel van de inverse functie £Î (als deze bestaat). 
Een vergelijking f(x) = a is ook te schrijven als £ Îa) z Xe 

In de tabel van de inverse functie is £Ìa) door interpolatie te 
benaderen, en zo krijgt men dus een benaderde oplossing van de 
gegeven vergelijking. Dit proces duidt men aan als inverse inter- 


polatie. 


Voor lineaire inverse interpolatie wordt de fout gegeven door: 


-hn 
ENE (y - vodty - va - h- met EE (yosyo + h). 


Probleem is dus het berekenen van de 2e afgeleide van de inverse 


functie. Ga, uitgaande van y = f(x), achtereenvolgens na 


d -1 EE 
ay If Ol = Fo) 
df ed gu 
dyz LE Cl == aen 3 


Bepaal aan de hand van onderstaande tabel m.b.v. lineaire inverse 


interpolatie het nulpunt van de functie In(x), bereken de gemaakte 


(Bn 


en controleer deze. 


„0. 10536 






0.18232 


Stel dat een voldoend vaak differentiëerbare functie f equidistant 


getabuleerd is op de punten 
XpsXg + h, Xx) + 2h, «… 


We vragen ons nu af hoe goed we nu de afgeleide van deze functie 
kunnen benaderen door alleen gebruik te maken van de getabuleerde 


waarden. 


a). Als eerste schatting voor de afgeleide van f in B Kot Ied 
bekijken we de eerste gedeelde differentie zoals die voorkomt 
in de volgende uitdrukking 


Ees) - fx) 


Ella, d = + R‚(k‚h) 





Toon aan dat voor de restterm geldt 


R‚(k‚h) = 5 £"(O) met x, SOS x 


k k+1 


b). Als tweede schatting voor de afgeleide van f in XX + k.h 
bezien we de eerste gedeelde differentie zoals die voorkomt in 


Eix dm EK ) 
OD) zt RO) > n 


Toon nu aan dat geldt 


n? 
R‚ hoh) = B £111(0) 


5 <0 << x 


“kel k+1 
ce). Teken eens een kromme lijn en beschouw deze als grafiek van een 
functie, Ga grafisch na welke van de 2 bovenstaande methoden 


betere resultaten doet verwachten 


d). 


e). 


E), 


7. 
We willen nu nagaan hoe dit in een concreet geval uitpakt, Stel, 
men heeft de sinus getabuleerd in de punten O‚,h,2h,...; dan 
kan men deze waarden gebruiken om de cosinus op b.v. [0,7/] 
te bepalen. 
Hoe groot mag de h nu nog zijn als U met de methode uit a) resp. 
b) de cosinus met een relatieve fout $ à. 1077 wilt berekenen? 
Uiteraard worden de fouten RK; & RK, kleiner als men h verkleint, 
en men verwacht derhalve zeer nauwkeurige uitkomsten als men h 
zeer klein neemt. Komt dit ook in de praktijk uit als U dit 


eens probeert op onderstaand tabelletje? 


x(rad) sin(x) cos{(x) 





Wat is de bron van de misère? 


Uit het voorafgaande concludeert men dat grote h schadelijk is 
wegens de R,‚ en R,, en dat kleine h schadelijk is wegens 

de afrondfouten in de tabelwaarden. We verwachten dus dat er 

een h is waarvoor de som van deze beide fouten minimaal zal zijn. 
Bepaal voor de methode onder a) en b) deze optimale h voorhet 
benaderen van de afgeleide van de sinus in het punt T/y als de 
sinus getabuleerd is in 3 cijfers achter de komma (zodat de 
onauwkeurigheid in de tabelwaarden hoogstens .….... is). U mag 
hierbij optredende factoren sin (0) en cos (0) wel door 


3/2 = 0.7 vervangen. Vergelijk de bereikbare nauwkeurigheden. 


8). 


Bij benaderingen speelt het begrip "orde van de benadering" een 
belangrijke rol, zoals we zagen : een proces met een foutterm 

z ch? levert voor h > 0 sneller nauwkeurige waarden op dan een 
proces met foutterm S c.h. We noemen de bedoelde processen resp. 
van de orde 2 en j. 

Het komt nogal eens voor dat men m.b.v. aldus bij benadering 
verkregen functiewaarden een nieuw benaderingsproces op gang zet, 
bijv. tussen deze waarden gaat interpoleren of gaat numeriek 
differentiëren. 

De vraag doet zich dan voor hoe de eerste benaderingsfouten in 
het tweede proces doorwerken. 

We illustreren dit m.b.v. de processen lineaire interpolatie en 


numerieke differentiatie. . 


Ga de volgende bewer ingen. zorgvuldig na 


- Wanneer men een C?-functie benadert door lineair interpoleren op 
x en xth dan heeft de benaderde functiewaarde een fout O(n?). 


= Wanneer we de afgeleide van f in x d.m.v. het proces 


tej Ex) » f(x+h) ope 


benaderen, vinden we f'(x) met een fout O(h?). 

- Wanneer bij lineaire interpolatie de functiewaarden met fouten 
Ee‚sE, belast zijn dan is het interpolatieresultaat hierdoor belast 
met een fout e,‚ zodat |el < max (le,l,le,l). 

… Wanneer bij numeriek differentiatie volgens procedé (* ) de functie- 
waarden met fouten EjsE, behept zijn, dan is het benaderde resul- 
taat hierdoor behept met een fout (e,-e‚)/2h. 

= Conclusie: Bij lineaire interpolatie zijn de fouten in de 
functiewaarden vrij onschuldig, bij numerieke differentiatie 


daarentegen worden de fouten voor kleinere h opgeblazen! 


II. 


9. 


Laat nu een tabel van functiewaarden gegeven zijn in de punten 

0, +h, +2h etc. 

We bepahen nu op 2 manieren een benadering voor de afgeleide in een 
niet-tabel punt x, gelegen tussen -jh en ih. 

Benader de afgeleiden in -Jh en jh m.b.v. procedé (*) uitgaande van 
£(O), flh) en f(-h). Deze worden benaderd met fout O(h?). 
Interpoleer nu lineair op de zojuist verkregen waarden om een be- 
nadering te verkrijgen voor de afgeleide in een punt x tussen -ih 
en +Èh. 


Het eindresultaat is dus behept met een fout O(h?). 


Uitgaande van 3 tabelpunten O‚,h en 2 h benaderen we f(x-th) door 
lineaire interpolatie op =-h en 0. Evenzo benaderen we f(xt}h) door 
lineaire interpolatie op O0 en h. Beide functiewaarden dus met een 
fout Oh), 

Bepaal nu m.b.v. GD) en fGed met een procedé analoog aan (*) 
een benadering voor f'(x). 

De totale fout in het eindresultaat is nu : O(h?) (vanwege de 


numerieke differentiatie) + O(h) (vanwege de onnauwkeurige functie- 


waarden) = O(h). 


Zoals we reeds in opgave 7 gezien hebben kan de eerste afgeleide 


in een punt x van een equidistante tabel behoorlijk benaderen door 
[f(xth) - flx-h)l/2h. 


We zullen nu een al eerder (nl. opgave 8) ontmoet procedé in 
concreto uitwerken. Stel we zoeken de afgeleide van een niet tabel- 
punt x, -jh < x < Jh, gebruikmakend van de functiewaarden f(-h), 


f(O) en £(h). 


10, 


e= 2 
a) Ga na : EM EKO gram) + B £trrE) 


He es 


b) Idem : ED ECR) = £'(-4h) + en LIINLE) 


ce) Interpoleer nu lineair tussen de linkerleden van a) en b). 


Ga na dat er komt 


xebh Flhj=fCOJ , Cphexh CECOjablsh) 


nh h h h 
2 
5 Ul Ef GerBh) Geej) + Be Etre) = 
z f'(x) + Olx). £'''(E) 
2 2 
met Dr Ss O(x) S =s s 


Wanneer worden deze grenzen aangenomen? 


10) In deze opgave willen we via enkele eenvoudige methoden de integraal 


1 
S e“ dx benaderen. 
0 


a). Benader de integraal met de trapeziumregel en vervolgens met 
de midpointregel. Bereken de fout t.o.v. het werkelijke antwoord. 
Wat valt U op bij vergelijking van het teken van de fout bij 
trapezium- resp. midpointregel? 
Was dit te verwachten? 

b). Pas ook eens de 1x gerepeteerde trapeziumregel en midpointregel 
toe. Vergelijk de verschillende fouten uit onderdeel a) en b) 
met elkaar en ga na of deze aan het door de theorie voorspelde 


gedrag voldoen. 





1.28403 
1.64872 
2.11700 
2.71828 








ec). Benader ook m.b.v. de Simpsonregel. 


ie pe 


12). 


Li 


vx 
Beschouw de integraal f S7 dx: 
1 





a). Ga na dat de integraal zeer snel naar 0 gaat voor Xx > w, 


Vermoeden is nu dat we de integraal redelijk kunen bena- 
M -x 





deren door f rs dx, M nader te bepalen. 
PR: —_ X 
b).Probeer nu f nr dx naar boven af te schatten. Omdat zowel 
je Er dx als 4 e * ax te bepalen zijn, zijn we nu geneigd 
a -M -X 
e e e 8 
rn door xr resp. door MZ te majoreren, en vervolgens te 
integreren. 


Van welke van de 2 methoden verwacht U de beste benadering? 
(doe dit op grond van het functiegedrag in ©). Toets Uw ver- 


moeden door expliciete berekening. 


ec). De procedure is nu duidelijk : Bij gegeven tolerantie e zullen 
 _-X 

we eerst M bepalen zodat : f Sr dx < 5 

M 


M %X 
en vervolgens een kwadratuurformule op { Sr dx toepassen 
1 


met tolerantie 6/9. 


Als men een integraal wil benaderen met de regel van Simpson 

kan men moeilijkheden verwachten als de vierde afgeleide van 

de integrand op het integratie interval oneindig wordt (waarom? ). 

Toch treedt in zo'n geval convergentie op als we de gerepeteerde 

Simpsonregel toepassen. We gaan na waarom dit zo is, en hoe snel 

die convergentie is. 

a). Toon convergentie aan, door de gerepeteerde Simpsonregel te 
interpreteren als Riemann sommen. 


D 


b). Leidt de regel van Simpson af voor f f(t)dt door een geschikte. 
a 
lineaire combinatie van de trapeziumregel en de midpoint- 


regel op [a,b]. Toon nu aan, dat de absolute waarde van de 





1e 


restterm van Simpson door de volgende uitdrukkingen ge- 


majoreerd wordt 





He) max _ |f'"(x)|, B max _ |f"'G)| 

1 xEl a ‚bl xEla,b) 
d 21 

ec). Hoeveel steunpunten zijn nodig om f x°*dx met de gerepeteerde 
0 


-h 
Simpsonregel te benaderen met een fout van hoogstens 10 7% 


En hoeveel steunpunten zijn minstens nodig als men de integraal 


benadert met een gerepeteerde midpointregel? 
d 


n / 
Resultaten verkregen bij werkelijke berekening van fx 2-dx met 
8 : À 


gerepeteerde Simpsonregel. 





n Xx gerept. Verhouding 





13). We bezien kwadratuurformules voor integralen van het type 
h 
S Yx. fh) dx 
0 


Zij fee A Dik. 


a). Ga aan de hand van de restterm na waarom de trapeziumregel 


het in het algemeen slecht zal doen. 


b). Nieuwe formule: Zij d,(x) het lineair interpolatiepolynoom 
van f met steunpunten O0 en h. 
h 
Integreer nu : f /x. d, Go) dx 
0 


ce). Laat zien dat de zo verkregen kwadratuurformule een fout 


: 7 
heeft van : 5 h 12 ECE Bx € (O‚h). 


d). Nog een formule: Zij b,(x) de som van de 0“ en 1° graadsterm 
h 


van de Taylorreeks van f rond 3/sh. Integreer : f /x d, Oeddx. 
0 
e). Bewijs dat de fout van de formule uit d) gelijk is aan 


HE em 
eh STE es ES (Oom. 





13. 


Ib. In deze opgave bezien we de volgende (singuliere) integraal 





1 e* 
I=f dx 
0 /x 
a) Ga na dat door de splitsing: 
1 EE: en 
af Sirel Sik 
0 /x 0 /x 


een elementaire integraal ontstaat en een integraal met een niet- 
singuliere integrand. Leent de laatste integraal zich voor toe- 
passen van de herhaalde trapeziumregel? En voor toepassing van 
de herhaalde midpointregel? 

b) Ga na dat we de singulariteit zonder problemen door kunnen 


drukken naar een hogere afgeleide b.v. 


1 da 
sf zi dx + f Bn dx. 
0 /x 0 /x 


(Het komt er dus op neer dat we een beginstuk van de Taylor- 
reeks van e” aftrekken). 
Dit zelfde kan ook door partiële integratie (ga na) 

ec) Een veel eleganteremethode is: substitueer t = /x. Ga na dat 
dan de integrand en z'n hogere afgeleiden geen singulariteiten 
bezitten. 
Hoe groot moet n zijn opdat dan een nx herhaalde trapeziumregel 


6 


een relatieve fout van hoogstens 10 ° oplevert? 


15. a) Laten we de berekende waarden van een functie f met absolute 


fouten $ e behept zijn. Laat zien dat dit het resultaat van 
b 


de (gerepeteerde) trapeziumregel toegepast op ff(x)dx belast met 
a 
een (extra) abs. fout S< (b-a) €. 


Laat zien dat dit bij de (gerepeteerde) midpoint en Simpson=- 


regels ook zo is. 


16. 


1u, 


b) Zie in dat dit het beste is wat te verwachten was, gezien de 
verandering van de ware integraal als men de functie f met € 


verstoort. 


In vraagstuk 15 zagen we dat abs. fouten $ e in de functiewaarden 
van een functie f aanleiding geven tot een (extra) abs. fout in 


1 
het gerepeteerde midpoint resultaat voor ff(x)dx van hoogstens €. 


In vraagstuk 14 pasten we de de so toe op 

() î El! dx. 

0 /x 

Stel dat de waarden van e”* opgeleverd worden met abs. fouten 
< e. Dan zijn de fouten in de functiewaarden die door de midpoint- 
regel gebruikt worden zeer veel groter dan €. 
We vragen ons af hoe de fouten in e“ doorwerken in het gerepeteerde 
midpoint resultaat. 
ad Zij I, het interval [(i-1)h, ihl, i = i,...‚n, met n = 1/,. 


Zie in dat de fout in e* een abs. fout in het midpointresultaat 


over Í, veroorzaakt $ £.v2h 
v21-1 
è Ei dx je 
b) Laat zien dat << Sf iP>2. 





VII=l iet (Zi 

ce) Laat nu zien dat de fout in e“ het gerepeteerde midpoint re- 
sultaat voor (* ) belast met een (extra) fout $ 2 €. 

d) Laat zien dat een fout van bijna 2e inderdaad kan optreden 
(zodat de foutschatting 2e wel realistisch genoemd mag worden); 
doe dit door de fouten in en geschikt te kiezen en een schatting 


analoog aan die in (b) echter met een >-teken te gebruiken. 








Geef een redelijke schatting voor h opdat Euler forward toegepast 
op Se = -2x een totale fout op het tijdstip t = 1.5 heeft die 


< 107* is. 


Gegeven de differentiaalvergelijking 


d% … 
dt = FEEN 


We beschouwen de volgende methode om deze vergelijking numeriek op 
te lossen: 
Xn 5 ax tm, ECE, ox )+B, ECE, 4 str). 
a) Laat zien dat er voor de differentiaalvergelijking: een = 0, 
X, = 1 geen convergentie optreedt als a # 1. 
b) Zij dus a = 1. Laat zien dat er voor de differentiaalvergelij- 
king = = 1, Xx = 0 geen convergentie optreedt als B,+B, # 1. 


ec) Uitgaande van de resultaten uit a) en b) luidt de methode in 


feite dus: 


Xn OF Kp-ath[BE,+C1-B)IE Jl. 


Bepaal de orde van deze methode als functie van de parameter B. 
d) Zij nu specifiek f(t‚x) = p.x(t)t+g(t), en de methode als in c). 
Bepaal de factor og zoals u die kent uit Hoofdstuk 3 8 2 e.v. 
e) Laat zien dat voor alle waarden van B geldt: 


gt/n 5 ePtsoqn) mits 0 SS tsS<T. 


We herinneren aan de volgende definitie: 


als er een constante K bestaat en een omgeving 


f = Olg) voor Xx > X, 


U van x, zodat: 
IÉGO)I S Klglx)l voor x E U. 
Het werken met O=symbolen is enerzijds handig, anderzijds moet men 


er wel mee uitkijken. 





a) 


b) 


ce) 


d) 


e) 


Er geldt voor vaste m en voor d = 1+ph dat 
GR = FPP 4 on) als h>0. 


Bewijs dit. 
Dit geldt echter niet uniform in m (d.w.z. de constante die in 
de O verscholen zit, stijgt met m boven alle grenzen) als p > 0. 


Laat dit zien. 


Laat zien dato" = elpn + O(h) wel geldt uniform in m als m 
voldoet aan mh $< T (T vast). ë 
d: 
Zij f € C'[o,T]. Laat zien dat hf(t,) = | f(s)ds + On?) 
Ted 
wiifern in Lals tx = ih @h Est. 5 [O.Tl, T vast. 
1 IR tE 

od ht 
Krachtens ec) geldt: h Z PCE. = | f(s)ds +O(h?), als 

' 1 

i=k+1 t 


k 


k en l vaste getallen zijn. 
Als echter het verschil tussen k en l toeneemt naarmate h af- 
neemt geldt dit niet meer. 


Toon aan dat: 


n tn 
h 5 Fltid = | f(s)ds + Oh) (*) 
ae 1 
1=0 0 
Ss 
als tn Te. 
Laat met behulp van de functie f(x) = e“ zien dat als men de 


eis tE, < T laat vallen zelfs de formule (*) uit d) niet meer 


behoeft te gelden. 


We beschouwen de differentiaalvergelijking 
OR oe e 
dt = 1000(1-x) 
x(0) = 1000 
N.B. Gebruik in deze opgave: e= 0.3679 


se men 


sm di 


a) 


b) 


e) 


d) 


e) 





Bepaal de oplossing. 
Bepaal m.b.v. Euler backward en Euler forward benaderde oplos- 
singswaarden voor: 


t = 0.001 met h = 0.001 


u 
oo 
Ld 
le} 
EN 


t = 0.01 met h 
t = 0.1 met h = 0.1 
welke formule geeft de beste resultaten? Had u dat al verwacht? 
Niettemin zijn ook de beste bij b) gevonden resultaten niet 
bijzonder goed, Dat zit hem daarin dat de oplossing zich in 
het begin vrijwel gedraagt als een oplossing van X' = -1000x. 
Een exacte oplossing hiervan wordt per stap h met an ver- 
menigvuldigd, terwijl onze methode dit geenszins doet. Wij 


trachten een methode te vinden die dat wel doet. We beschouwen 


daartoe de methode: 


X s OK # h[BFCt, ‚x,)+(1-B)f(t )] 


n+1 n n+1°*n+1 


en proberen B zo te bepalen dat de benaderde oplossing van 


-1000h 


xXx! = =1000x ook met een factor e per stap wordt vermenig- 


vuldigd. Ga na dat volgt: 


e”1000h(444000n)-1 


B SAL, 
-1000h(1-e" PUR, 


Los de oorspronkelijke differentiaalvergelijking nu ook nog 


eens benaderd op m.b.v. de nieuwe methode en wel 


in t = 0.001 met h = 0.001 (u dient te vinden 8 = 0.419), 
en in t = 0.01 met h = 0.01 (laat zien dat 8 = 5 en re- 
ken dan verder met 8 = ee 


De moraal van dit verhaal is dat Euler backward wel aardig is 
als we eenmaal bij die waarden van t zijn beland waar de op- 
lossing langzaam varieert, maar dat deze methode ontoereikend 


is om in het begin een snel variërend gedrag goed op te leveren 


met een toch nog wat royale stap. 
De methode met een aan dit snel variërende gedrag aangepaste 8 


(men spreekt van "exponential fitting") doet het dan veel beter. 


TE eeen 


1. Xen onderzoekt vaak de convergentie van een successief-substi ee 
proces aan de hand van een prent, waarbij men al op andere gronden 
weet dat in een zeker interval precies één nulpunt van de functie 
LieË, 

Kies als vergelijking: x?°+2x+1 = 0. 

a) Laat in de eerste plaats zien dat deze vergelijking een enkel- 
voudige wortel heeft op (=150). 

b) We schrijven de vergelijking nu op 2 verschillende wijzen als 


Xx = lx) en wel met 


2xt1 1 
Cx) = dx) = TDEXE 


Teken nu voor b, en d, een successief substitutie plaatje. 
Controleer of u de plaatjes enigszins redelijk getekend hebt met 
name in die gebieden waar het van belang is dat KEE Se Ar 

e) Ga na bij welke van de twee processen er sprake is van conver- 
gentie naar een vast punt, en geef dan een boven- en ondergrens 
voor de convergentiefactor die hoogstens een factor 2 verschil- 


len. 


2. We bezien de vergelijking e* = atyx; a>1, 0S<y<t1. 
a) Laat zien dat er precies één positieve wortel a 18, 
b) Zie eerst grafisch en daarna ook analytisch in dat het itera- 


tieve proces 
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voor geen positive startwaarde (# a) convergeert naar a. 


ce) Laat zien dat het proces 


x loglatyx,) Ie els Bs an 


k+1 
naar a convergeert voor elke positieve startwaarde x,. 


d) Laat zien dat het proces uit ec) lineair convergent is met 





(asymptotische) convergentiefactor le. 
e) Laat zien dat de rij der iteranden in ec) monotoon dalend res- 
pectievelijk monotoon stijgend is voor x, > a resp. Xg S a. 


f) Ga na dat geldt: 
log(a) < a < log(2a). 


Hint: gebruik achtereenvolgens als weartiaanden 0 en log(2a) 
en maak gebruik van log(2a) < a (hetgeen u ook nog even bewij- 
zen moet). 

g) Laat zien dat 


y/2a < convergentiefactor < y/a. 


Stel men wenst te bepalen x = E, voor gegeven a (a f 0) zonder dat 
men de mogelijkheid heeft een deling "echt" uit te voeren (in de 
begintijd van de rekenmachines was dit van belang). 


a) Ga na dat N.R. toegepast op a-d = 0 het volgende proces oplevert: 


Xx ee 2x;-axs Ls Bass Cx) 


wed 
en hierin wordt niet meer gedeeld. 

b) Ga zowel grafisch als aan de hand van de recursie uit a) na 
dat er convergentie optreedt voor 0 < x, S 2, 


e) Zie in dat er monotone convergentie optreedt als 0 < x, S 5 
wat gebeurt er als À Hr 22 
d) Wat ons nu nog rest om het proces uit a) toe te kunnen passen 


is een "redelijke" startwaarde X9. 


d dan kan 1074 als start- 


Een mogelijkheid is: als 109"? <a < 10 
waarde gebruikt worden. 
e) Ga na dat a asymptotisch als factor tussen hij en hi optreedt. 
f) De eigenschap uit e betekent echter niet dat, om een zekere re- 


latieve nauwkeurigheid te behalen, men in het geval van grote 


a veel meer iteratieslagen nodig zou hebben dan in het geval 














van kleine a. 
Ga hiertoe m.b.v. e) na dat de relatieve fout bij iedere slag 
ongeveer gekwadrateerd wordt. 
Conelusie: voor alle waarden van a gaat het even goed. 

g) De conclusie uit f) was ook al te halen uit formule (+) uit a), 
vermenigvuldigd met a aan beide kanten. 
Zie in dat hieruit volgt dat de relatieve fout in Xj uitslui- 


tend afhangt van i en de relatieve fout in x,. 
a) Toon aan dat voor 


EG) == xPax2ad la Sdj. 


het Newton-Raphson proces met startwaarde x, = 10° aanvankelijk 


lineair convergeert met als factor © nd 





b) Laat zien dat voor 
F(x) = e*-b b EeR 


het Newton-Raphson proces, voor voldoend grote xs, willekeurig 
langzaam kan convergeren, d.w.z. dat hs 44/h; willekeurig dicht 


bij 1 kan liggen. 


We willen m.b.v. N.R. een tabel van een inverse functie gaan maken. 
we kiezen hiervoor arcsin, die we getabelleerd willen hebber voor 


xE [0,3/2] met stapjes h. De op te lossen vgl. luidt dus 
sin(y)=jh = 0. e (+) 
Voor 1 = 0 is deze triviaal oplosbaar. 
a) Zij a de wortel van (+). Laat zien dat N.R. gegarandeerd con- 
vergeert naar a; voor (+) wanneer men als startwaarde Fa ° ed 


kiest. 


B} Zij h; = Vids bij het in a) bedoelde proces. Laat zien dat 


[hs‚4/hil <3. 


i+1 





c) 


We 


a) 


b) 


Stel nu dat men van het in a) bedoelde proces voor iedere j 
slechts 1 stap uitvoert, zodat men dus de waarde van y, als 
benaderde waarde voor arcsin(ih) aflevert. 

Laat nu met volledig inductie zien dat de fout in de aldus 

opgeleverde waarde van de arcsin kleiner dan h? is, als nog 


gegeven is h $< 0.1. 


willen aantonen dat een iteratief proces met orde kleiner dan 


geen convergent proces kan zijn. 


Laat eerst zien: de: k; = ch; en [ks 44! > Ek 


(Zie bewijs Stelling 3.9.) 


Zie nu in dat er inderdaad geen convergentie kan optreden. 


ie. 
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2e DEELTENTAMEN TOEGEPASTE ANALYSE 


ZATERDAG 14 DECEMBER 10,00 — 13.00 UUR 


- Zet op elk blad uw naam. 





-— Zet op het dubbele vel, dat verder niet beschreven wordt, en 

| waarin het werk ingeleverd wordt, uw naam en adres en studie- 
richting. 

- Maak, zo vaak u het maar nodig acht, gebruik van het aanwij=- 


zingenblad. 





EE 





Leid uitgaande van een interpolatiepolynoom af dat de fout van de 


trapeziumregel toegepast op B flxhde, ES C*la,bl,is 


as 
n (ba) £"'(E) ES Db. 


Stel dat men cos(x) tabelleert met stapgrootte h op het interval 


Lon] . Hoe groot moet men h dan nemen opdat de fout bij lineaire 
interpolatie $ Ee is? Zelfde vraag voor kwadratische interpo- 
latie. 


Laat zien dat het n x gerepeteerde midpoint-resultaat toegepast 
op Ek fFGddx, f E C[a,bl, naar de ware waarde van de integraal 
convergeert voor n *@, (N.B. Als dat u gemakkelijker is mag u 


zelfs wel aannemen f € C* [a‚b}). 


Bepaal een redelijke schatting voor het onzekerheidsinterval voor 


het nulpunt van f(x) = e“-3, als f(x) met een abs. fout S dap 
berekend wordt. 
Beschouw de vergelijking e“-3 = 0. Laat m.b.v. een plaatje zien 


dat het Newton-Raphson-proces toegepast op e*-3 voor elke start- 
waarde naar de wortel convergeert. Deel iets mede over monotonie 
eigenschappen van de opgeleverde rij iteranden (in afhankelijk- 
heid van de startwaarde). Is er sprake van kwadratische conver- 


gentie? 


a) Bepaal voor de kwadratuurformule: 
ofP Fax = wjfCO) + wzfCh) + waf'(0) + wf!) + RON) 
de gewichten Wi zó dat polynomen van zo hoog mogelijke graad 


nog exact door deze formule geïntegreerd worden. 


b) Aan welke formule doet de gevonden kwadratuurformule u denken? 
c) Bepaal c en k in de formule voor R(h), aannemend dat deze de 
gedaante heeft 
A nd CP (0 <E <n) 
d) Toegift: Kunt u aantonen dat de restterm de in ec) genoemde ge- 


daante heeft? 


Geef aan de hand van de (gerepeteerde) Simpson-regel een korte en 


duidelijk uiteenzetting over "glijdende stap". 
8 g 








Aanwijzingenblad 


Zij gegeven: Xjs--:»X, uit het interval [a,b 


Posse sP, niet-negatieve gehele getallen 


een functie f die minstens max P; maal differentieer- 
| 1 
ĳ baar is. 


Als het interpolatiepolynoom p(x) van f(x) voldoet aan: 


(k) = (k) e 
P (xj) z= f xj) voor k = Oj.s.,P; 


LS Uri 


(N+1) 


dan geldt als N = ntryt...tr,, en als f bestaat op (a,b) dat 


er voor alle x € [a,bl een £ € (a,b) bestaat zodat 
(x-X es iP 
En n N+1 . 


roti ri+1 


Elx) = px) == xx) Cx-xj) 


Voor de midpoint-regel. geldt de volgende relatie 


3 
SSP ear = flats) + Ip £"E), E € (a,ath) 


als f € C*[a,athl. 








